A Counterfactual Simulation Model of Causal Judgment

1 Setup

1.1 Load packages

library("knitr")     # for knitting stuff
library("lme4")      # for linear mixed effects models
library("broom")     # for tidying model results 
library("brms")      # Bayesian regression with Stan
library("corrr")     # for tidy correlation matrix
library("xtable")    # for latex tables
library("jsonlite")  # for reading json files
library("janitor")   # cleans stuff
library("ggtern")    # for ternary plots
library("Hmisc")     # bootstrapped confidence intervals
library("tidybayes") # for Bayesian data analysis
library("png")       # adding pngs to images
library("grid")      # functions for dealing with images 
library("egg")       # for geom_custom
library("tidyverse") # for wrangling, plotting, etc. 

1.2 Helper functions

# root mean squared error
rmse = function(x, y){
  return(sqrt(mean((x - y)^2)))
}

actual_counterfactual_plot = function(x, ylabel) {
  p = ggplot(data = x, 
              mapping = aes(x = clipindex,
                            y = value,
                            fill = colorindex,
                            group = model)) +
    stat_summary(geom = "bar",
                 fun.y = "mean",
                 color = "black",
                 position = position_dodge(0.7),
                 width = 0.7) +
    stat_summary(geom = "linerange",
                 fun.data = "mean_cl_boot",
                 position = position_dodge(0.7),
                 size = 1) +
    geom_point(data = x %>% mutate(value = ifelse(model == "rating", value, NA)),
               position = position_jitterdodge(dodge.width = 0.7, jitter.width = 0.2),
               shape = 21,
               color = "black",
               size = 1,
               alpha = 0.2) +
    scale_fill_manual(values = c("red2", "green2", "black")) +
    facet_grid(index_actual ~ index_counterfactual) +
    geom_text(data = df.text %>% filter(!is.na(label)),
              mapping = aes(x = x, y = y, label = as.character(label)),
              size = 6,
              position = position_dodge(width = .7),
              hjust = 0.5) +
    coord_cartesian(xlim = c(0.5, 2.5),
                    ylim = c(-5, 105),
                    expand = F,
                    clip = "off") +
    labs(y = ylabel) +
    theme_bw() +
    theme(
      text = element_text(size = 20),
      legend.position = "none",
      axis.text.x = element_blank(),
      axis.ticks.x = element_blank(),
      axis.title.y = element_text(margin = margin(0, 10, 0, 0)),
      panel.spacing.y = unit(c(.8), "cm"),
      plot.margin = unit(c(0.2, 0.2, .8, .2), "cm"),
      axis.title.x = element_blank(),
      panel.grid = element_blank(),
      strip.background = element_blank(),
      panel.border = element_rect(colour = "black")
    )
  print(p)
}

actual_counterfactual_threeballs_plot = function(x, ylabel) {
  p = ggplot(data = x, 
              mapping = aes(x = ball,
                            y = value,
                            group = model,
                            fill = colorindex)) +
    stat_summary(fun.y = "mean",
                 geom = "bar",
                 color = "black",
                 position = position_dodge(0.9),
                 width = 0.9) +
    stat_summary(fun.data = "mean_cl_boot",
                 geom = "linerange",
                 color = "black",
                 position = position_dodge(0.9)) +
    geom_point(data = x %>% mutate(value = ifelse(model == "rating", value, NA)),
               position = position_jitterdodge(jitter.width = 0.3, jitter.height = 0, dodge.width = 0.9),
               shape = 21,
               color = "black",
               size = 1,
               alpha = 0.15) +
    facet_wrap(~clip, ncol = 8) +
    geom_text(data = df.text, 
              mapping = aes(x = x,
                            y = y,
                            label = as.character(label)),
              size = 5, 
              position = position_dodge(width = .9), 
              hjust = 0.5, 
              fontface = "bold") +
    coord_cartesian(xlim = c(0.5, 2.5), 
                    ylim = c(-5, 120), 
                    expand = F,
                    clip = "off") +
    scale_y_continuous(breaks = seq(0, 100, 25)) +
    scale_fill_manual(values = c("red2", "green2", "black")) +
    scale_color_manual(values = c("red2", "green2", "black")) +
    labs(y = ylabel) +
    theme_bw() +
    theme(
      text = element_text(size = 16),
      legend.position = "none",
      axis.title.x = element_blank(),
      panel.spacing.y = unit(c(.3), "cm"),
      panel.grid = element_blank(),
      strip.background = element_blank(),
      strip.text = element_blank(),
      panel.border = element_rect(colour = "black")
    )
  print(p)
}

# set ggplot theme 
theme_set(
  theme_classic()
)

2 Experiment 1: 2 balls

2.1 Clips

include_graphics("../../figures/diagrams/exp1_diagrams.png")

2.2 Read in data

# causal judgments
df.exp1.causal = read.csv("../../data/experiment1_causal.csv", 
                          header = T,
                          stringsAsFactors = F)

# counterfactual judgments
df.exp1.counterfactual = read.csv("../../data/experiment1_counterfactual.csv",
                                  header = T,
                                  stringsAsFactors = F)

# Information about each clip 
df.exp1.clipinfo = tibble(
  clip = 1:18,
  outcome_actual = c(0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1),
  outcome_counterfactual = c(0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1),
  index_actual = rep(c("actual miss", 
                       "actual close", 
                       "actual hit"), each = 6),
  index_counterfactual = rep(rep(c("counterfactual miss",
                                   "counterfactual close", 
                                   "counterfactual hit"), each = 2), 3)
) %>% 
  mutate(index_actual = factor(index_actual, levels = c("actual miss", 
                                                        "actual close", 
                                                        "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss",
                                                  "counterfactual close", 
                                                  "counterfactual hit")))

# Structure data 

df.exp1.causal.long = df.exp1.causal$data %>%
  strsplit("\n") %>%
  unlist() %>%
  strsplit(",") %>%
  unlist() %>%
  matrix(ncol = 36, byrow = T) %>%
  as.data.frame(stringsAsFactors = F) %>%
  mutate_all("as.numeric") %>%
  setNames(str_c(c("clip_", "rating_"), rep(1:18, each = 2))) %>%
  mutate(participant = 1:nrow(.)) %>%
  gather("index", "value", -participant) %>% 
  separate(index, into = c("index", "clip")) %>% 
  spread(index, value) %>% 
  left_join(df.exp1.clipinfo,
            by = "clip") %>%
  arrange(participant, clip) %>%
  ungroup()

df.exp1.counterfactual.long = df.exp1.counterfactual$data %>%
  strsplit("\n") %>%
  unlist() %>%
  strsplit(",") %>%
  unlist() %>%
  matrix(ncol = 36, byrow = T) %>%
  as.data.frame(stringsAsFactors = F) %>%
  mutate_all("as.numeric") %>%
  setNames(str_c(c("clip_", "rating_"), rep(1:18, each = 2))) %>%
  mutate(participant = 1:nrow(.)) %>%
  gather("index", "value", -participant) %>% 
  separate(index, into = c("index", "clip")) %>% 
  spread(index, value) %>% 
  left_join(df.exp1.clipinfo,
            by = "clip") %>%
  arrange(participant, clip) %>%
  ungroup()

2.2.1 Demographics

# counterfactual condition 
df.exp1.counterfactual %>% 
  separate(extra,
           into = c("gender", "age", "difficulty", "smoothness", "time"),
           sep = ",") %>% 
  mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>% 
  select(gender:time) %>% 
  mutate_if(is.character, ~ as.numeric(.)) %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2))

# causal condition 
df.exp1.causal %>% 
  separate(extra,
           into = c("gender", "age", "difficulty", "smoothness", "time"),
           sep = ",") %>% 
  mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>% 
  select(gender:time) %>% 
  mutate_if(is.character, ~ as.numeric(.)) %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2))

   n age_mean age_sd n_female time_mean time_sd
1 41       32    8.8       19       7.5    3.12
   n age_mean age_sd n_female time_mean time_sd
1 41       33   10.3       21      9.65    6.58

2.3 Read in model predictions

2.3.1 Counterfactual condition

# read model predictions 
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "2ball")]

df.exp1.model = files %>%
  map(~ read.csv(str_c(path, .))) %>%
  reduce(rbind) %>%
  set_names(c("clip", "noise", "prediction")) %>%
  left_join(df.exp1.clipinfo, by = "clip") %>%
  mutate(prediction = ifelse(outcome_actual == 1, 1 - prediction, prediction))

# calculate mean counterfactual judgments 
df.exp1.counterfactual.means = df.exp1.counterfactual.long %>%
  group_by(clip) %>%
  summarize(rating_mean = mean(rating),
            rating_low = smean.cl.boot(rating)[2],
            rating_high = smean.cl.boot(rating)[3]) %>%
  ungroup() %>%
  left_join(df.exp1.clipinfo, by = "clip")

# find noisy simulation model that correlates best with the mean counterfactual judgments
df.exp1.counterfactual.model = df.exp1.model %>% 
  group_by(noise) %>% 
  nest() %>% 
  mutate(correlation = map(data, ~ cor(.x$prediction,
                                       df.exp1.counterfactual.means$rating_mean)),
         correlation = unlist(correlation)) %>% 
  ungroup() %>% 
  filter(correlation == max(correlation)) %>% 
  unnest(data) %>% 
  mutate(prediction = prediction * 100)

2.3.2 Causal condition

df.exp1.causal.model = df.exp1.counterfactual.long %>% 
  group_by(clip) %>% 
  summarize(empirical = mean(rating)) %>% 
  ungroup() %>% 
  left_join(df.exp1.counterfactual.model %>% 
              select(-c(noise, correlation)),
            by = "clip") %>% 
  rename(simulation = prediction) %>% 
  mutate(empirical = ifelse(outcome_actual == 1, 100 - empirical, empirical),
         simulation = ifelse(outcome_actual == 1, 100 - simulation, simulation))

2.4 Counterfactual condition

2.4.1 Plots

df.plot = df.exp1.counterfactual.long %>% 
  mutate(rating = abs(rating),
         clipindex = rep(c(1, 2), nrow(.)/2)) %>% 
  left_join(df.exp1.counterfactual.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual")) %>% 
  gather("model", "value", c(rating, prediction)) %>% 
  mutate(colorindex = ifelse(outcome_counterfactual == 1, 2, 1),
         colorindex = ifelse(model != "rating", 3, colorindex),
         colorindex = as.factor(colorindex),
         index_actual = factor(index_actual,levels = c("actual miss", "actual close", "actual hit")),
         index_counterfactual = factor(index_counterfactual,levels = c("counterfactual miss", "counterfactual close", "counterfactual hit")),
         model = factor(model, levels = c("rating", "prediction"))) %>% 
  mutate(clipindex = ifelse(clip == 11, 2, clipindex),
         clipindex = ifelse(clip == 12, 1, clipindex)) #swap clips 11 and 12


df.text = df.plot %>%
  expand(index_actual, index_counterfactual, clipindex, model) %>%
  mutate(
    label = rep(1:18, each = 2),
    label = ifelse(model != "rating", NA, label),
    y = -15,
    x = clipindex,
    colorindex = NA
  )

actual_counterfactual_plot(df.plot, "counterfactual judgment")

Warning: Removed 738 rows containing missing values (geom_point).

# ggsave("../../figures/plots/exp1_counterfactual_bars.pdf",
#        width = 8,
#        height = 6)

2.4.2 Stats

2.4.2.1 Model evaluation

df.exp1.counterfactual.means %>% 
  left_join(df.exp1.counterfactual.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual")) %>% 
  summarize(simulation_r = cor(rating_mean, prediction),
            simulation_rmse = rmse(rating_mean, prediction)) %>% 
  mutate_all(~ round(., 2))

# A tibble: 1 x 2
  simulation_r simulation_rmse
         <dbl>           <dbl>
1         0.98            12.0

2.5 Causal condition

2.5.1 Plots

model.name = c("empirical")
# model.name = c("simulation")

df.plot = df.exp1.causal.long %>%
  mutate(
    rating = abs(rating),
    clipindex = rep(c(1, 2), nrow(.) / 2)
  ) %>%
  left_join(df.exp1.causal.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual")) %>%
  gather("model", "value", c(rating, simulation, empirical)) %>%
  filter(model %in% c(model.name, "rating")) %>%
  mutate(model = factor(model, levels = c("rating", model.name))) %>%
  mutate(
    colorindex = ifelse(outcome_actual == 1, 2, 1),
    colorindex = ifelse(model != "rating", 3, colorindex),
    colorindex = as.factor(colorindex),
    index.actual = factor(index_actual, 
                          levels = c("actual miss", "actual close", "actual hit")),
    index.counterfactual = factor(index_counterfactual, 
                                  levels = c("counterfactual miss", 
                                             "counterfactual close", 
                                             "counterfactual hit"))
  ) %>%
  mutate(
    clipindex = ifelse(clip == 11, 2, clipindex),
    clipindex = ifelse(clip == 12, 1, clipindex)
  ) # swap clips 11 and 12

df.text = df.plot %>%
  expand(index_actual, index_counterfactual, clipindex, model) %>%
  mutate(
    label = rep(1:18, each = 2),
    label = ifelse(model != "rating", NA, label),
    y = -15,
    x = clipindex,
    colorindex = NA
  )

actual_counterfactual_plot(df.plot, "causal judgment")

Warning: Removed 738 rows containing missing values (geom_point).

# ggsave("../../figures/plots/exp1_causal_bars.pdf",
#        width = 8,
#        height = 6)

2.5.2 Stats

2.5.2.1 Model evaluation

df.exp1.causal.long %>% 
  mutate(rating = abs(rating)) %>% 
  group_by(clip,
           outcome_actual, 
           outcome_counterfactual,
           index_actual,
           index_counterfactual) %>% 
  summarize(rating = mean(rating)) %>% 
  ungroup() %>% 
  left_join(df.exp1.causal.model,
            by = c("clip", 
                   "outcome_actual", 
                   "outcome_counterfactual", 
                   "index_actual", 
                   "index_counterfactual")) %>% 
  summarize(empirical_r = cor(rating, empirical),
            empirical_rmse = rmse(rating, empirical),
            simulation_r = cor(rating, simulation),
            simulation_rmse = rmse(rating, simulation)) %>% 
  mutate_all(~ round(., 2))

# A tibble: 1 x 4
  empirical_r empirical_rmse simulation_r simulation_rmse
        <dbl>          <dbl>        <dbl>           <dbl>
1        0.96           8.57         0.94            16.8

2.5.2.2 Bayesian regression

# random intercepts model 
fit.brm_exp1_causal = brm(formula = rating ~ 1 + index_actual + index_counterfactual + outcome_actual + (1 + index_actual + index_counterfactual + outcome_actual | participant),
                          data = df.exp1.causal.long,
                          cores = 2,
                          seed = 1,
                          file = "cache/brm_exp1_causal")

fit.brm_exp1_causal %>% 
  summary()

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: rating ~ 1 + index_actual + index_counterfactual + outcome_actual + (1 + index_actual + index_counterfactual + outcome_actual | participant) 
   Data: df.exp1.causal.long (Number of observations: 738) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~participant (Number of levels: 41) 
                                                                                   Estimate
sd(Intercept)                                                                         13.95
sd(index_actualactualclose)                                                            3.40
sd(index_actualactualhit)                                                              5.51
sd(index_counterfactualcounterfactualclose)                                           11.02
sd(index_counterfactualcounterfactualhit)                                             21.28
sd(outcome_actual)                                                                    22.11
cor(Intercept,index_actualactualclose)                                                -0.18
cor(Intercept,index_actualactualhit)                                                  -0.03
cor(index_actualactualclose,index_actualactualhit)                                     0.11
cor(Intercept,index_counterfactualcounterfactualclose)                                 0.15
cor(index_actualactualclose,index_counterfactualcounterfactualclose)                   0.03
cor(index_actualactualhit,index_counterfactualcounterfactualclose)                     0.03
cor(Intercept,index_counterfactualcounterfactualhit)                                  -0.43
cor(index_actualactualclose,index_counterfactualcounterfactualhit)                     0.10
cor(index_actualactualhit,index_counterfactualcounterfactualhit)                       0.17
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit)    -0.05
cor(Intercept,outcome_actual)                                                         -0.48
cor(index_actualactualclose,outcome_actual)                                           -0.02
cor(index_actualactualhit,outcome_actual)                                             -0.26
cor(index_counterfactualcounterfactualclose,outcome_actual)                           -0.17
cor(index_counterfactualcounterfactualhit,outcome_actual)                             -0.44
                                                                                   Est.Error
sd(Intercept)                                                                           3.12
sd(index_actualactualclose)                                                             2.59
sd(index_actualactualhit)                                                               4.10
sd(index_counterfactualcounterfactualclose)                                             4.46
sd(index_counterfactualcounterfactualhit)                                               3.91
sd(outcome_actual)                                                                      3.93
cor(Intercept,index_actualactualclose)                                                  0.37
cor(Intercept,index_actualactualhit)                                                    0.36
cor(index_actualactualclose,index_actualactualhit)                                      0.38
cor(Intercept,index_counterfactualcounterfactualclose)                                  0.29
cor(index_actualactualclose,index_counterfactualcounterfactualclose)                    0.37
cor(index_actualactualhit,index_counterfactualcounterfactualclose)                      0.36
cor(Intercept,index_counterfactualcounterfactualhit)                                    0.19
cor(index_actualactualclose,index_counterfactualcounterfactualhit)                      0.36
cor(index_actualactualhit,index_counterfactualcounterfactualhit)                        0.34
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit)      0.28
cor(Intercept,outcome_actual)                                                           0.18
cor(index_actualactualclose,outcome_actual)                                             0.36
cor(index_actualactualhit,outcome_actual)                                               0.36
cor(index_counterfactualcounterfactualclose,outcome_actual)                             0.28
cor(index_counterfactualcounterfactualhit,outcome_actual)                               0.19
                                                                                   l-95% CI
sd(Intercept)                                                                          8.14
sd(index_actualactualclose)                                                            0.14
sd(index_actualactualhit)                                                              0.18
sd(index_counterfactualcounterfactualclose)                                            2.13
sd(index_counterfactualcounterfactualhit)                                             14.23
sd(outcome_actual)                                                                    15.17
cor(Intercept,index_actualactualclose)                                                -0.79
cor(Intercept,index_actualactualhit)                                                  -0.70
cor(index_actualactualclose,index_actualactualhit)                                    -0.64
cor(Intercept,index_counterfactualcounterfactualclose)                                -0.41
cor(index_actualactualclose,index_counterfactualcounterfactualclose)                  -0.70
cor(index_actualactualhit,index_counterfactualcounterfactualclose)                    -0.67
cor(Intercept,index_counterfactualcounterfactualhit)                                  -0.75
cor(index_actualactualclose,index_counterfactualcounterfactualhit)                    -0.63
cor(index_actualactualhit,index_counterfactualcounterfactualhit)                      -0.52
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit)    -0.60
cor(Intercept,outcome_actual)                                                         -0.78
cor(index_actualactualclose,outcome_actual)                                           -0.68
cor(index_actualactualhit,outcome_actual)                                             -0.84
cor(index_counterfactualcounterfactualclose,outcome_actual)                           -0.66
cor(index_counterfactualcounterfactualhit,outcome_actual)                             -0.76
                                                                                   u-95% CI
sd(Intercept)                                                                         20.41
sd(index_actualactualclose)                                                            9.69
sd(index_actualactualhit)                                                             15.38
sd(index_counterfactualcounterfactualclose)                                           19.72
sd(index_counterfactualcounterfactualhit)                                             29.22
sd(outcome_actual)                                                                    30.76
cor(Intercept,index_actualactualclose)                                                 0.59
cor(Intercept,index_actualactualhit)                                                   0.65
cor(index_actualactualclose,index_actualactualhit)                                     0.78
cor(Intercept,index_counterfactualcounterfactualclose)                                 0.70
cor(index_actualactualclose,index_counterfactualcounterfactualclose)                   0.70
cor(index_actualactualhit,index_counterfactualcounterfactualclose)                     0.70
cor(Intercept,index_counterfactualcounterfactualhit)                                  -0.02
cor(index_actualactualclose,index_counterfactualcounterfactualhit)                     0.73
cor(index_actualactualhit,index_counterfactualcounterfactualhit)                       0.77
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit)     0.47
cor(Intercept,outcome_actual)                                                         -0.08
cor(index_actualactualclose,outcome_actual)                                            0.68
cor(index_actualactualhit,outcome_actual)                                              0.50
cor(index_counterfactualcounterfactualclose,outcome_actual)                            0.41
cor(index_counterfactualcounterfactualhit,outcome_actual)                             -0.06
                                                                                   Rhat
sd(Intercept)                                                                      1.00
sd(index_actualactualclose)                                                        1.00
sd(index_actualactualhit)                                                          1.00
sd(index_counterfactualcounterfactualclose)                                        1.00
sd(index_counterfactualcounterfactualhit)                                          1.00
sd(outcome_actual)                                                                 1.01
cor(Intercept,index_actualactualclose)                                             1.00
cor(Intercept,index_actualactualhit)                                               1.00
cor(index_actualactualclose,index_actualactualhit)                                 1.00
cor(Intercept,index_counterfactualcounterfactualclose)                             1.00
cor(index_actualactualclose,index_counterfactualcounterfactualclose)               1.00
cor(index_actualactualhit,index_counterfactualcounterfactualclose)                 1.00
cor(Intercept,index_counterfactualcounterfactualhit)                               1.00
cor(index_actualactualclose,index_counterfactualcounterfactualhit)                 1.01
cor(index_actualactualhit,index_counterfactualcounterfactualhit)                   1.01
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) 1.00
cor(Intercept,outcome_actual)                                                      1.01
cor(index_actualactualclose,outcome_actual)                                        1.02
cor(index_actualactualhit,outcome_actual)                                          1.01
cor(index_counterfactualcounterfactualclose,outcome_actual)                        1.00
cor(index_counterfactualcounterfactualhit,outcome_actual)                          1.00
                                                                                   Bulk_ESS
sd(Intercept)                                                                          1599
sd(index_actualactualclose)                                                            1553
sd(index_actualactualhit)                                                              1127
sd(index_counterfactualcounterfactualclose)                                             805
sd(index_counterfactualcounterfactualhit)                                              1162
sd(outcome_actual)                                                                     1065
cor(Intercept,index_actualactualclose)                                                 3783
cor(Intercept,index_actualactualhit)                                                   4114
cor(index_actualactualclose,index_actualactualhit)                                     1953
cor(Intercept,index_counterfactualcounterfactualclose)                                 1567
cor(index_actualactualclose,index_counterfactualcounterfactualclose)                    828
cor(index_actualactualhit,index_counterfactualcounterfactualclose)                     1345
cor(Intercept,index_counterfactualcounterfactualhit)                                   1034
cor(index_actualactualclose,index_counterfactualcounterfactualhit)                      379
cor(index_actualactualhit,index_counterfactualcounterfactualhit)                        497
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit)      740
cor(Intercept,outcome_actual)                                                           805
cor(index_actualactualclose,outcome_actual)                                             360
cor(index_actualactualhit,outcome_actual)                                               432
cor(index_counterfactualcounterfactualclose,outcome_actual)                             900
cor(index_counterfactualcounterfactualhit,outcome_actual)                              1742
                                                                                   Tail_ESS
sd(Intercept)                                                                          1756
sd(index_actualactualclose)                                                            1582
sd(index_actualactualhit)                                                              1691
sd(index_counterfactualcounterfactualclose)                                            1162
sd(index_counterfactualcounterfactualhit)                                              2131
sd(outcome_actual)                                                                     1768
cor(Intercept,index_actualactualclose)                                                 3088
cor(Intercept,index_actualactualhit)                                                   2881
cor(index_actualactualclose,index_actualactualhit)                                     2573
cor(Intercept,index_counterfactualcounterfactualclose)                                 2576
cor(index_actualactualclose,index_counterfactualcounterfactualclose)                   2032
cor(index_actualactualhit,index_counterfactualcounterfactualclose)                     2298
cor(Intercept,index_counterfactualcounterfactualhit)                                   2121
cor(index_actualactualclose,index_counterfactualcounterfactualhit)                      897
cor(index_actualactualhit,index_counterfactualcounterfactualhit)                       1619
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit)     1765
cor(Intercept,outcome_actual)                                                          1632
cor(index_actualactualclose,outcome_actual)                                             810
cor(index_actualactualhit,outcome_actual)                                              1213
cor(index_counterfactualcounterfactualclose,outcome_actual)                            1574
cor(index_counterfactualcounterfactualhit,outcome_actual)                              2685

Population-Level Effects: 
                                        Estimate Est.Error l-95% CI
Intercept                                 -18.33      3.36   -24.86
index_actualactualclose                     4.31      3.40    -2.54
index_actualactualhit                       4.06      4.84    -5.51
index_counterfactualcounterfactualclose   -17.52      3.37   -24.36
index_counterfactualcounterfactualhit     -61.37      4.38   -69.95
outcome_actual                             92.36      5.11    82.21
                                        u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept                                 -11.66 1.00     2026     2429
index_actualactualclose                    10.94 1.00     2989     2747
index_actualactualhit                      13.27 1.00     2531     2706
index_counterfactualcounterfactualclose   -10.84 1.00     2898     2310
index_counterfactualcounterfactualhit     -52.89 1.00     2210     2724
outcome_actual                            102.38 1.00     2179     2885

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    30.09      0.89    28.42    31.83 1.00     3208     3142

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

2.5.2.3 Specific hypotheses

df.exp1.causal.posterior = df.exp1.causal.long %>% 
  group_by(clip, index_actual, index_counterfactual, outcome_actual) %>% 
  summarize(value = mean(rating)) %>% 
  ungroup() %>% 
  add_fitted_draws(newdata = .,
                   model = fit.brm_exp1_causal,
                   re_formula = NA) %>% 
  ungroup() %>% 
  select(clip, .value, .draw) %>% 
  spread(clip, .value)

func_posterior_difference = function(data, clip1, clip2){
  data %>% 
    mutate(difference = .data[[clip1]] - .data[[clip2]]) %>% 
    pull(difference) %>% 
    mean_hdci() %>% 
    mutate_if(is.numeric, ~round(., 2)) %>% 
    summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])")) %>% 
    print()  
}

func_posterior_difference(data = df.exp1.causal.posterior,
                          clip1 = "8",
                          clip2 = "13")

            difference
1 (0.26 [-6.51, 7.31])

2.5.3 Tables

fit.brm_exp1_causal %>% 
  tidy() %>% 
  filter(str_detect(term, "b_")) %>% 
  mutate(term = str_remove(term, "b_"),
         term = tolower(term)) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  rename(`lower 95% HDI` = lower,
         `upper 95% HDI` = upper) %>% 
  xtable() %>%
  print(include.rownames = F,
        booktabs = T)

% latex table generated in R 3.6.1 by xtable 1.8-4 package
% Tue Nov 12 15:05:19 2019
\begin{table}[ht]
\centering
\begin{tabular}{lrrrr}
  \toprule
term & estimate & std.error & lower 95\% HDI & upper 95\% HDI \\ 
  \midrule
intercept & -18.33 & 3.36 & -23.81 & -12.78 \\ 
  index\_actualactualclose & 4.31 & 3.40 & -1.40 & 9.83 \\ 
  index\_actualactualhit & 4.06 & 4.84 & -4.01 & 11.88 \\ 
  index\_counterfactualcounterfactualclose & -17.52 & 3.37 & -22.96 & -12.01 \\ 
  index\_counterfactualcounterfactualhit & -61.37 & 4.38 & -68.54 & -54.30 \\ 
  outcome\_actual & 92.36 & 5.11 & 84.04 & 100.96 \\ 
   \bottomrule
\end{tabular}
\end{table}

3 Experiment 2: Brick & Teleport

3.1 Clips

include_graphics("../../figures/diagrams/exp2_diagrams.png")

3.2 Read in data

df.exp2.causal_first = read.csv("../../data/experiment2_causal_first.csv", header = T, stringsAsFactors = F)
df.exp2.counterfactual_first = read.csv("../../data/experiment2_counterfactual_first.csv", header = T, stringsAsFactors = F)

# Information about each clip

df.exp2.clipinfo = tibble(
    clip = 1:20,
    outcome_actual = c(0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1),
    outcome_counterfactual = c(0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1),
    index_actual = c(rep(c("actual miss", "actual close", "actual hit"), each = 6), rep("practice", 2)),
    index_counterfactual = c(rep(rep(c("counterfactual miss", "counterfactual close", "counterfactual hit"), each = 2), 3), rep("practice", 2))
  ) %>% 
  mutate(index_actual = factor(index_actual, levels = c("actual miss", 
                                                        "actual close", 
                                                        "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss",
                                                  "counterfactual close", 
                                                  "counterfactual hit")))

# Structure data
df.exp2.causal_first.long = df.exp2.causal_first$data %>%
  strsplit("\n") %>%
  unlist() %>%
  strsplit(",") %>%
  unlist() %>%
  matrix(ncol = 80, byrow = T) %>%
  as.data.frame(stringsAsFactors = F) %>%
  mutate_all("as.numeric") %>%
  setNames(str_c(c("clip_", "rating_"), rep(1:40, each = 2))) %>%
  mutate(participant = 1:nrow(.)) %>%
  gather("index", "value", -participant) %>% 
  separate(index, into = c("index", "trial")) %>% 
  mutate(trial = as.numeric(trial)) %>% 
  spread(index, value) %>% 
  arrange(participant, trial) %>%
  mutate(
    question = rep(rep(c("causal", "counterfactual"), each = 20), max(participant)),
    condition = "causal_first"
  ) %>%
  left_join(df.exp2.clipinfo,
            by = "clip") %>%
  arrange(participant, question, clip) %>%
  select(participant, question, clip, rating, everything())

df.exp2.counterfactual_first.long = df.exp2.counterfactual_first$data %>%
  strsplit("\n") %>%
  unlist() %>%
  strsplit(",") %>%
  unlist() %>%
  matrix(ncol = 80, byrow = T) %>%
  as.data.frame(stringsAsFactors = F) %>%
  mutate_all("as.numeric") %>%
  setNames(str_c(c("clip_", "rating_"), rep(1:40, each = 2))) %>%
  mutate(participant = 1:nrow(.)) %>%
  gather("index", "value", -participant) %>% 
  separate(index, into = c("index", "trial")) %>% 
  mutate(trial = as.numeric(trial)) %>% 
  spread(index, value) %>% 
  arrange(participant, trial) %>%
  mutate(
    question = rep(rep(c("counterfactual", "causal"), each = 20), max(participant)),
    condition = "counterfactual_first"
  ) %>%
  left_join(df.exp2.clipinfo,
            by = "clip") %>% 
  arrange(participant, question, clip) %>%
  select(participant, question, clip, rating, everything())

# combine data from both conditions
df.exp2.combined.long = df.exp2.causal_first.long %>%
  bind_rows(df.exp2.counterfactual_first.long %>% 
              mutate(participant = participant + max(df.exp2.causal_first.long$participant)))

3.2.1 Demographics

# counterfactual condition 
df.exp2.counterfactual_first %>% 
  separate(extra,
           into = c("gender", "age", "difficulty", "smoothness", "time", "condition"),
           sep = ",") %>% 
  select(gender:time) %>% 
  bind_rows(df.exp2.causal_first %>% 
  separate(extra,
           into = c("gender", "age", "difficulty", "smoothness", "time", "condition"),
           sep = ",") %>% 
  select(gender:time)) %>% 
  mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>% 
  mutate_if(is.character, ~ as.numeric(.)) %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2))

   n age_mean age_sd n_female time_mean time_sd
1 82       35   12.2       45     21.25    5.11

3.3 Read in model predictions

3.3.1 Counterfactual condition

# read model predictions 
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "teleport")]

df.exp2.model = files %>%
  map(~ read.csv(str_c(path, .))) %>% 
  reduce(rbind) %>%
  set_names(c("clip", "noise", "prediction")) %>%
  left_join(df.exp2.clipinfo, by = "clip") %>%
  mutate(prediction = ifelse(outcome_actual == 1, 1 - prediction, prediction))

# calculate mean counterfactual judgments 
df.exp2.counterfactual.means = df.exp2.combined.long %>%
  filter(question == "counterfactual",
         clip <= 18) %>%
  # mutate(rating = abs(rating)) %>% 
  group_by(clip) %>%
  summarize(rating = mean(rating)) %>%
  ungroup() %>%
  left_join(df.exp2.clipinfo, by = "clip")

# find noisy simulation model that correlates best with the mean counterfactual judgments
df.exp2.counterfactual.model = df.exp2.model %>% 
  group_by(noise) %>% 
  nest() %>% 
  mutate(correlation = map(data, ~ cor(.x$prediction,
                                       df.exp2.counterfactual.means$rating)),
         correlation = unlist(correlation)) %>% 
  ungroup() %>% 
  filter(correlation == max(correlation)) %>% 
  unnest(data) %>% 
  mutate(prediction = prediction * 100)

3.3.2 Causal condition

# Model predictions based on counterfactual judgments
df.exp2.causal.model = df.exp2.combined.long %>%
  filter(
    question == "counterfactual",
    clip <= 18
  ) %>%
  group_by(clip, condition) %>%
  summarize(rating = mean(rating)) %>%
  spread(condition, rating) %>% 
  left_join(df.exp2.combined.long %>%
              filter(
                question == "counterfactual",
                clip <= 18
              ) %>%
              group_by(clip) %>%
              summarize(combined = mean(rating)),
            by = "clip") %>% 
  left_join(df.exp2.counterfactual.model,
            by = "clip") %>% 
  select(-c(correlation, noise)) %>% 
  rename(simulation = prediction) %>% 
  mutate_at(.vars = vars(simulation, causal_first, counterfactual_first, combined), .funs = list(~ ifelse(outcome_actual == 1, 100 - ., .)))

3.4 Counterfactual condition

3.4.1 Plots

df.plot = df.exp2.combined.long %>%
  filter(question == "counterfactual",
         clip <= 18) %>%
  mutate(clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
  left_join(df.exp2.counterfactual.model,
            by = c("clip", 
                   "outcome_actual", 
                   "outcome_counterfactual", 
                   "index_actual", 
                   "index_counterfactual")) %>%
  gather("model", "value", c(rating, prediction)) %>%
  mutate(colorindex = ifelse(outcome_counterfactual == 1, 2, 1),
         colorindex = ifelse(model != "rating", 3, colorindex),
         colorindex = as.factor(colorindex),
         index_actual = factor(index_actual,
                               levels = c("actual miss", 
                                          "actual close", 
                                          "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss", 
                                                  "counterfactual close", 
                                                  "counterfactual hit")),
         model = factor(model,
                             levels = c("rating", "prediction")))

df.text = df.plot %>%
  expand(index_actual, index_counterfactual, clipindex, model) %>%
  mutate(
    label = rep(c(
      "1 (A)", "2 (B)", "3", "4", "5 (A)", "6 (B)",
      "7", "8 (C)", "9", "10", "11", "12 (C)",
      "13 (D)", "14 (E)", "15", "16", "17 (D)", "18 (E)"
    ), each = 2),
    label = ifelse(model != "rating", NA, label),
    y = -15,
    x = clipindex,
    colorindex = NA
  )

actual_counterfactual_plot(df.plot, "counterfactual judgment")

# ggsave("../../figures/plots/exp2_counterfactual_bars.pdf",
#        width = 8,
#        height = 6)

3.4.2 Stats

3.4.2.1 Model evaluation

df.exp2.counterfactual.means %>% 
  left_join(df.exp2.counterfactual.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual")) %>% 
  summarize(simulation_r = cor(rating, prediction),
            simulation_rmse = rmse(rating, prediction)) %>% 
  mutate_all(~ round(., 2))

# A tibble: 1 x 2
  simulation_r simulation_rmse
         <dbl>           <dbl>
1         0.97            14.6

3.4.2.2 Bayesian regressions

3.4.2.2.1 Check for potential order effects

# random intercepts & slopes model 
fit.brm_exp2_counterfactual = brm(formula = rating ~ 1 + condition * (index_actual + index_counterfactual) + (1 + index_actual + index_counterfactual | participant),
                                  data = df.exp2.combined.long %>% 
                                    na.omit() %>% 
                                    filter(question == "counterfactual"),
                                  cores = 2,
                                  seed = 1,
                                  file = "cache/brm_exp2_counterfactual")

fit.brm_exp2_counterfactual %>% 
  summary()

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: rating ~ 1 + condition * (index_actual + index_counterfactual) + (1 + index_actual + index_counterfactual | participant) 
   Data: df.exp2.combined.long %>% na.omit() %>% filter(que (Number of observations: 1476) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~participant (Number of levels: 82) 
                                                                                   Estimate
sd(Intercept)                                                                          8.80
sd(index_actualactualclose)                                                            6.76
sd(index_actualactualhit)                                                              2.37
sd(index_counterfactualcounterfactualclose)                                           14.65
sd(index_counterfactualcounterfactualhit)                                             22.12
cor(Intercept,index_actualactualclose)                                                 0.13
cor(Intercept,index_actualactualhit)                                                  -0.02
cor(index_actualactualclose,index_actualactualhit)                                     0.18
cor(Intercept,index_counterfactualcounterfactualclose)                                -0.41
cor(index_actualactualclose,index_counterfactualcounterfactualclose)                  -0.04
cor(index_actualactualhit,index_counterfactualcounterfactualclose)                    -0.05
cor(Intercept,index_counterfactualcounterfactualhit)                                  -0.86
cor(index_actualactualclose,index_counterfactualcounterfactualhit)                    -0.20
cor(index_actualactualhit,index_counterfactualcounterfactualhit)                      -0.11
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit)     0.36
                                                                                   Est.Error
sd(Intercept)                                                                           1.84
sd(index_actualactualclose)                                                             3.07
sd(index_actualactualhit)                                                               1.78
sd(index_counterfactualcounterfactualclose)                                             2.43
sd(index_counterfactualcounterfactualhit)                                               2.57
cor(Intercept,index_actualactualclose)                                                  0.31
cor(Intercept,index_actualactualhit)                                                    0.40
cor(index_actualactualclose,index_actualactualhit)                                      0.42
cor(Intercept,index_counterfactualcounterfactualclose)                                  0.20
cor(index_actualactualclose,index_counterfactualcounterfactualclose)                    0.29
cor(index_actualactualhit,index_counterfactualcounterfactualclose)                      0.38
cor(Intercept,index_counterfactualcounterfactualhit)                                    0.08
cor(index_actualactualclose,index_counterfactualcounterfactualhit)                      0.27
cor(index_actualactualhit,index_counterfactualcounterfactualhit)                        0.39
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit)      0.16
                                                                                   l-95% CI
sd(Intercept)                                                                          5.13
sd(index_actualactualclose)                                                            0.58
sd(index_actualactualhit)                                                              0.10
sd(index_counterfactualcounterfactualclose)                                            9.96
sd(index_counterfactualcounterfactualhit)                                             17.37
cor(Intercept,index_actualactualclose)                                                -0.46
cor(Intercept,index_actualactualhit)                                                  -0.74
cor(index_actualactualclose,index_actualactualhit)                                    -0.69
cor(Intercept,index_counterfactualcounterfactualclose)                                -0.73
cor(index_actualactualclose,index_counterfactualcounterfactualclose)                  -0.61
cor(index_actualactualhit,index_counterfactualcounterfactualclose)                    -0.73
cor(Intercept,index_counterfactualcounterfactualhit)                                  -0.97
cor(index_actualactualclose,index_counterfactualcounterfactualhit)                    -0.70
cor(index_actualactualhit,index_counterfactualcounterfactualhit)                      -0.78
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit)     0.01
                                                                                   u-95% CI
sd(Intercept)                                                                         12.58
sd(index_actualactualclose)                                                           12.43
sd(index_actualactualhit)                                                              6.67
sd(index_counterfactualcounterfactualclose)                                           19.56
sd(index_counterfactualcounterfactualhit)                                             27.51
cor(Intercept,index_actualactualclose)                                                 0.74
cor(Intercept,index_actualactualhit)                                                   0.74
cor(index_actualactualclose,index_actualactualhit)                                     0.85
cor(Intercept,index_counterfactualcounterfactualclose)                                 0.02
cor(index_actualactualclose,index_counterfactualcounterfactualclose)                   0.55
cor(index_actualactualhit,index_counterfactualcounterfactualclose)                     0.71
cor(Intercept,index_counterfactualcounterfactualhit)                                  -0.66
cor(index_actualactualclose,index_counterfactualcounterfactualhit)                     0.37
cor(index_actualactualhit,index_counterfactualcounterfactualhit)                       0.66
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit)     0.63
                                                                                   Rhat
sd(Intercept)                                                                      1.00
sd(index_actualactualclose)                                                        1.01
sd(index_actualactualhit)                                                          1.00
sd(index_counterfactualcounterfactualclose)                                        1.00
sd(index_counterfactualcounterfactualhit)                                          1.00
cor(Intercept,index_actualactualclose)                                             1.00
cor(Intercept,index_actualactualhit)                                               1.00
cor(index_actualactualclose,index_actualactualhit)                                 1.00
cor(Intercept,index_counterfactualcounterfactualclose)                             1.00
cor(index_actualactualclose,index_counterfactualcounterfactualclose)               1.00
cor(index_actualactualhit,index_counterfactualcounterfactualclose)                 1.01
cor(Intercept,index_counterfactualcounterfactualhit)                               1.00
cor(index_actualactualclose,index_counterfactualcounterfactualhit)                 1.00
cor(index_actualactualhit,index_counterfactualcounterfactualhit)                   1.00
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) 1.00
                                                                                   Bulk_ESS
sd(Intercept)                                                                          1720
sd(index_actualactualclose)                                                             516
sd(index_actualactualhit)                                                              1448
sd(index_counterfactualcounterfactualclose)                                            1135
sd(index_counterfactualcounterfactualhit)                                              1621
cor(Intercept,index_actualactualclose)                                                 1418
cor(Intercept,index_actualactualhit)                                                   3693
cor(index_actualactualclose,index_actualactualhit)                                     2473
cor(Intercept,index_counterfactualcounterfactualclose)                                  887
cor(index_actualactualclose,index_counterfactualcounterfactualclose)                    474
cor(index_actualactualhit,index_counterfactualcounterfactualclose)                      331
cor(Intercept,index_counterfactualcounterfactualhit)                                    652
cor(index_actualactualclose,index_counterfactualcounterfactualhit)                     1422
cor(index_actualactualhit,index_counterfactualcounterfactualhit)                       1323
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit)     2216
                                                                                   Tail_ESS
sd(Intercept)                                                                          1686
sd(index_actualactualclose)                                                             639
sd(index_actualactualhit)                                                              2198
sd(index_counterfactualcounterfactualclose)                                            1936
sd(index_counterfactualcounterfactualhit)                                              2274
cor(Intercept,index_actualactualclose)                                                 2112
cor(Intercept,index_actualactualhit)                                                   3102
cor(index_actualactualclose,index_actualactualhit)                                     2861
cor(Intercept,index_counterfactualcounterfactualclose)                                 1737
cor(index_actualactualclose,index_counterfactualcounterfactualclose)                    771
cor(index_actualactualhit,index_counterfactualcounterfactualclose)                      714
cor(Intercept,index_counterfactualcounterfactualhit)                                   1299
cor(index_actualactualclose,index_counterfactualcounterfactualhit)                     1430
cor(index_actualactualhit,index_counterfactualcounterfactualhit)                       2815
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit)     3157

Population-Level Effects: 
                                                                      Estimate
Intercept                                                                13.38
conditioncounterfactual_first                                            -1.11
index_actualactualclose                                                  -1.48
index_actualactualhit                                                    -4.42
index_counterfactualcounterfactualclose                                  37.94
index_counterfactualcounterfactualhit                                    73.80
conditioncounterfactual_first:index_actualactualclose                    -0.55
conditioncounterfactual_first:index_actualactualhit                      -0.85
conditioncounterfactual_first:index_counterfactualcounterfactualclose     1.23
conditioncounterfactual_first:index_counterfactualcounterfactualhit      -0.02
                                                                      Est.Error
Intercept                                                                  2.58
conditioncounterfactual_first                                              3.65
index_actualactualclose                                                    2.78
index_actualactualhit                                                      2.40
index_counterfactualcounterfactualclose                                    3.33
index_counterfactualcounterfactualhit                                      4.17
conditioncounterfactual_first:index_actualactualclose                      3.84
conditioncounterfactual_first:index_actualactualhit                        3.41
conditioncounterfactual_first:index_counterfactualcounterfactualclose      4.79
conditioncounterfactual_first:index_counterfactualcounterfactualhit        5.92
                                                                      l-95% CI
Intercept                                                                 8.30
conditioncounterfactual_first                                            -7.99
index_actualactualclose                                                  -6.79
index_actualactualhit                                                    -9.15
index_counterfactualcounterfactualclose                                  31.43
index_counterfactualcounterfactualhit                                    65.54
conditioncounterfactual_first:index_actualactualclose                    -8.29
conditioncounterfactual_first:index_actualactualhit                      -7.60
conditioncounterfactual_first:index_counterfactualcounterfactualclose    -8.37
conditioncounterfactual_first:index_counterfactualcounterfactualhit     -11.96
                                                                      u-95% CI
Intercept                                                                18.38
conditioncounterfactual_first                                             6.04
index_actualactualclose                                                   3.95
index_actualactualhit                                                     0.20
index_counterfactualcounterfactualclose                                  44.60
index_counterfactualcounterfactualhit                                    81.99
conditioncounterfactual_first:index_actualactualclose                     6.79
conditioncounterfactual_first:index_actualactualhit                       5.87
conditioncounterfactual_first:index_counterfactualcounterfactualclose    10.75
conditioncounterfactual_first:index_counterfactualcounterfactualhit      11.59
                                                                      Rhat
Intercept                                                             1.00
conditioncounterfactual_first                                         1.00
index_actualactualclose                                               1.00
index_actualactualhit                                                 1.00
index_counterfactualcounterfactualclose                               1.00
index_counterfactualcounterfactualhit                                 1.00
conditioncounterfactual_first:index_actualactualclose                 1.00
conditioncounterfactual_first:index_actualactualhit                   1.00
conditioncounterfactual_first:index_counterfactualcounterfactualclose 1.00
conditioncounterfactual_first:index_counterfactualcounterfactualhit   1.00
                                                                      Bulk_ESS
Intercept                                                                 2268
conditioncounterfactual_first                                             1966
index_actualactualclose                                                   2898
index_actualactualhit                                                     3619
index_counterfactualcounterfactualclose                                   2118
index_counterfactualcounterfactualhit                                     1818
conditioncounterfactual_first:index_actualactualclose                     2827
conditioncounterfactual_first:index_actualactualhit                       3265
conditioncounterfactual_first:index_counterfactualcounterfactualclose     2006
conditioncounterfactual_first:index_counterfactualcounterfactualhit       1684
                                                                      Tail_ESS
Intercept                                                                 2863
conditioncounterfactual_first                                             2564
index_actualactualclose                                                   2790
index_actualactualhit                                                     3303
index_counterfactualcounterfactualclose                                   2627
index_counterfactualcounterfactualhit                                     2490
conditioncounterfactual_first:index_actualactualclose                     3028
conditioncounterfactual_first:index_actualactualhit                       3250
conditioncounterfactual_first:index_counterfactualcounterfactualclose     2617
conditioncounterfactual_first:index_counterfactualcounterfactualhit       2692

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    26.57      0.56    25.50    27.72 1.00     2894     2705

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

3.4.3 Tables

fit.brm_exp2_counterfactual %>% 
  tidy() %>% 
  filter(str_detect(term, "b_")) %>% 
  mutate(term = str_remove(term, "b_"),
         term = tolower(term)) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  rename(`lower 95% HDI` = lower,
         `upper 95% HDI` = upper) %>% 
  xtable() %>%
  print(include.rownames = F,
        booktabs = T)

% latex table generated in R 3.6.1 by xtable 1.8-4 package
% Tue Nov 12 15:05:24 2019
\begin{table}[ht]
\centering
\begin{tabular}{lrrrr}
  \toprule
term & estimate & std.error & lower 95\% HDI & upper 95\% HDI \\ 
  \midrule
intercept & 13.38 & 2.58 & 9.12 & 17.55 \\ 
  conditioncounterfactual\_first & -1.11 & 3.65 & -6.99 & 4.82 \\ 
  index\_actualactualclose & -1.48 & 2.78 & -6.02 & 3.11 \\ 
  index\_actualactualhit & -4.42 & 2.40 & -8.32 & -0.46 \\ 
  index\_counterfactualcounterfactualclose & 37.94 & 3.33 & 32.39 & 43.35 \\ 
  index\_counterfactualcounterfactualhit & 73.80 & 4.17 & 67.00 & 80.57 \\ 
  conditioncounterfactual\_first:index\_actualactualclose & -0.55 & 3.84 & -7.11 & 5.73 \\ 
  conditioncounterfactual\_first:index\_actualactualhit & -0.85 & 3.41 & -6.52 & 4.85 \\ 
  conditioncounterfactual\_first:index\_counterfactualcounterfactualclose & 1.23 & 4.79 & -6.62 & 9.18 \\ 
  conditioncounterfactual\_first:index\_counterfactualcounterfactualhit & -0.02 & 5.92 & -9.86 & 9.55 \\ 
   \bottomrule
\end{tabular}
\end{table}

3.5 Causal condition

3.5.1 Plots

# which causal condition? 
# condition.name = c("causal_first")
condition.name = c("counterfactual_first")
# condition.name = c("causal_first", "counterfactual_first")

# which model for the counterfactual simulations 
# model.name = c("simulation")
# model.name = c("causal_first")
# model.name = c("counterfactual_first")
model.name = c("combined")

df.plot = df.exp2.combined.long %>%
  filter(
    question == "causal",
    clip <= 18
  ) %>%
  filter(condition %in% condition.name) %>%
  mutate(
    rating = abs(rating),
    clipindex = rep(c(1, 2), nrow(.) / 2)
  ) %>%
  left_join(df.exp2.causal.model,
            by = c("clip", 
                   "outcome_actual", 
                   "outcome_counterfactual", 
                   "index_actual", 
                   "index_counterfactual")) %>%
  gather("model", "value", c("rating", model.name)) %>%
  mutate(model = factor(model, levels = c("rating", model.name))) %>%
  mutate(
    colorindex = ifelse(outcome_actual == 1, 2, 1),
    colorindex = ifelse(model != "rating", 3, colorindex),
    colorindex = as.factor(colorindex),
    index_actual = factor(index_actual, 
                          levels = c("actual miss", 
                                     "actual close", 
                                     "actual hit")),
    index_counterfactual = factor(index_counterfactual, 
                                  levels = c("counterfactual miss", 
                                             "counterfactual close", 
                                             "counterfactual hit")))

df.text = df.plot %>%
  expand(index_actual, index_counterfactual, clipindex, model) %>%
  mutate(
    label = rep(c(
      "1 (A)", "2 (B)", "3", "4", "5 (A)", "6 (B)",
      "7", "8 (C)", "9", "10", "11", "12 (C)",
      "13 (D)", "14 (E)", "15", "16", "17 (D)", "18 (E)"
    ), each = 2),
    label = ifelse(model != "rating", NA, label),
    y = -15,
    x = clipindex,
    colorindex = NA
  )

actual_counterfactual_plot(df.plot, "causal judgment")

Warning: Removed 738 rows containing missing values (geom_point).

# ggsave(str_c("../../figures/plots/exp2_", condition.name ,"_causal_bars.pdf"),
#        width = 8,
#        height = 6)

3.5.2 Stats

3.5.2.1 Model evaluation

df.exp2.combined.long %>% 
  filter(question == "causal",
         clip <= 18) %>% 
  mutate(rating = abs(rating)) %>% 
  group_by(clip,
           outcome_actual, 
           outcome_counterfactual,
           index_actual,
           index_counterfactual) %>% 
  summarize(rating = mean(rating)) %>% 
  ungroup() %>% 
  left_join(df.exp2.causal.model,
            by = c("clip", 
                   "outcome_actual", 
                   "outcome_counterfactual", 
                   "index_actual", 
                   "index_counterfactual")) %>% 
  summarize(empirical_r = cor(rating, combined),
            empirical_rmse = rmse(rating, combined),
            simulation_r = cor(rating, simulation),
            simulation_rmse = rmse(rating, simulation)) %>% 
  mutate_all(~ round(., 2))

# A tibble: 1 x 4
  empirical_r empirical_rmse simulation_r simulation_rmse
        <dbl>          <dbl>        <dbl>           <dbl>
1        0.97           8.36         0.94            19.4

3.5.2.2 Bayesian regression

# random intercepts + slopes model 
fit.brm_exp2_causal = brm(formula = rating ~ 1 + condition * (index_actual + index_counterfactual + outcome_actual) +
                            (1 + index_actual + index_counterfactual + outcome_actual | participant),
                          data = df.exp2.combined.long %>% 
                            filter(question == "causal",
                                   clip <= 18),
                          cores = 2,
                          seed = 1,
                          file = "cache/brm_exp2_causal")

fit.brm_exp2_causal %>% 
  summary()

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: rating ~ 1 + condition * (index_actual + index_counterfactual + outcome_actual) + (1 + index_actual + index_counterfactual + outcome_actual | participant) 
   Data: df.exp2.combined.long %>% filter(question == "caus (Number of observations: 1476) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~participant (Number of levels: 82) 
                                                                                   Estimate
sd(Intercept)                                                                         15.12
sd(index_actualactualclose)                                                            3.57
sd(index_actualactualhit)                                                              6.72
sd(index_counterfactualcounterfactualclose)                                           15.91
sd(index_counterfactualcounterfactualhit)                                             22.85
sd(outcome_actual)                                                                    23.85
cor(Intercept,index_actualactualclose)                                                -0.16
cor(Intercept,index_actualactualhit)                                                  -0.43
cor(index_actualactualclose,index_actualactualhit)                                     0.18
cor(Intercept,index_counterfactualcounterfactualclose)                                 0.02
cor(index_actualactualclose,index_counterfactualcounterfactualclose)                   0.19
cor(index_actualactualhit,index_counterfactualcounterfactualclose)                     0.03
cor(Intercept,index_counterfactualcounterfactualhit)                                  -0.62
cor(index_actualactualclose,index_counterfactualcounterfactualhit)                     0.18
cor(index_actualactualhit,index_counterfactualcounterfactualhit)                       0.39
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit)     0.20
cor(Intercept,outcome_actual)                                                         -0.69
cor(index_actualactualclose,outcome_actual)                                            0.09
cor(index_actualactualhit,outcome_actual)                                              0.02
cor(index_counterfactualcounterfactualclose,outcome_actual)                           -0.13
cor(index_counterfactualcounterfactualhit,outcome_actual)                              0.07
                                                                                   Est.Error
sd(Intercept)                                                                           2.20
sd(index_actualactualclose)                                                             2.41
sd(index_actualactualhit)                                                               3.79
sd(index_counterfactualcounterfactualclose)                                             2.75
sd(index_counterfactualcounterfactualhit)                                               2.88
sd(outcome_actual)                                                                      2.96
cor(Intercept,index_actualactualclose)                                                  0.35
cor(Intercept,index_actualactualhit)                                                    0.31
cor(index_actualactualclose,index_actualactualhit)                                      0.37
cor(Intercept,index_counterfactualcounterfactualclose)                                  0.19
cor(index_actualactualclose,index_counterfactualcounterfactualclose)                    0.34
cor(index_actualactualhit,index_counterfactualcounterfactualclose)                      0.32
cor(Intercept,index_counterfactualcounterfactualhit)                                    0.11
cor(index_actualactualclose,index_counterfactualcounterfactualhit)                      0.35
cor(index_actualactualhit,index_counterfactualcounterfactualhit)                        0.31
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit)      0.18
cor(Intercept,outcome_actual)                                                           0.12
cor(index_actualactualclose,outcome_actual)                                             0.34
cor(index_actualactualhit,outcome_actual)                                               0.30
cor(index_counterfactualcounterfactualclose,outcome_actual)                             0.18
cor(index_counterfactualcounterfactualhit,outcome_actual)                               0.17
                                                                                   l-95% CI
sd(Intercept)                                                                         11.02
sd(index_actualactualclose)                                                            0.16
sd(index_actualactualhit)                                                              0.48
sd(index_counterfactualcounterfactualclose)                                           10.82
sd(index_counterfactualcounterfactualhit)                                             17.57
sd(outcome_actual)                                                                    18.35
cor(Intercept,index_actualactualclose)                                                -0.75
cor(Intercept,index_actualactualhit)                                                  -0.87
cor(index_actualactualclose,index_actualactualhit)                                    -0.59
cor(Intercept,index_counterfactualcounterfactualclose)                                -0.34
cor(index_actualactualclose,index_counterfactualcounterfactualclose)                  -0.53
cor(index_actualactualhit,index_counterfactualcounterfactualclose)                    -0.61
cor(Intercept,index_counterfactualcounterfactualhit)                                  -0.80
cor(index_actualactualclose,index_counterfactualcounterfactualhit)                    -0.56
cor(index_actualactualhit,index_counterfactualcounterfactualhit)                      -0.36
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit)    -0.17
cor(Intercept,outcome_actual)                                                         -0.88
cor(index_actualactualclose,outcome_actual)                                           -0.60
cor(index_actualactualhit,outcome_actual)                                             -0.56
cor(index_counterfactualcounterfactualclose,outcome_actual)                           -0.47
cor(index_counterfactualcounterfactualhit,outcome_actual)                             -0.26
                                                                                   u-95% CI
sd(Intercept)                                                                         19.73
sd(index_actualactualclose)                                                            8.90
sd(index_actualactualhit)                                                             14.42
sd(index_counterfactualcounterfactualclose)                                           21.60
sd(index_counterfactualcounterfactualhit)                                             28.81
sd(outcome_actual)                                                                    30.03
cor(Intercept,index_actualactualclose)                                                 0.56
cor(Intercept,index_actualactualhit)                                                   0.32
cor(index_actualactualclose,index_actualactualhit)                                     0.80
cor(Intercept,index_counterfactualcounterfactualclose)                                 0.41
cor(index_actualactualclose,index_counterfactualcounterfactualclose)                   0.78
cor(index_actualactualhit,index_counterfactualcounterfactualclose)                     0.64
cor(Intercept,index_counterfactualcounterfactualhit)                                  -0.37
cor(index_actualactualclose,index_counterfactualcounterfactualhit)                     0.78
cor(index_actualactualhit,index_counterfactualcounterfactualhit)                       0.86
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit)     0.52
cor(Intercept,outcome_actual)                                                         -0.42
cor(index_actualactualclose,outcome_actual)                                            0.71
cor(index_actualactualhit,outcome_actual)                                              0.61
cor(index_counterfactualcounterfactualclose,outcome_actual)                            0.23
cor(index_counterfactualcounterfactualhit,outcome_actual)                              0.38
                                                                                   Rhat
sd(Intercept)                                                                      1.00
sd(index_actualactualclose)                                                        1.00
sd(index_actualactualhit)                                                          1.00
sd(index_counterfactualcounterfactualclose)                                        1.00
sd(index_counterfactualcounterfactualhit)                                          1.00
sd(outcome_actual)                                                                 1.00
cor(Intercept,index_actualactualclose)                                             1.00
cor(Intercept,index_actualactualhit)                                               1.00
cor(index_actualactualclose,index_actualactualhit)                                 1.00
cor(Intercept,index_counterfactualcounterfactualclose)                             1.00
cor(index_actualactualclose,index_counterfactualcounterfactualclose)               1.01
cor(index_actualactualhit,index_counterfactualcounterfactualclose)                 1.01
cor(Intercept,index_counterfactualcounterfactualhit)                               1.00
cor(index_actualactualclose,index_counterfactualcounterfactualhit)                 1.01
cor(index_actualactualhit,index_counterfactualcounterfactualhit)                   1.02
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) 1.00
cor(Intercept,outcome_actual)                                                      1.00
cor(index_actualactualclose,outcome_actual)                                        1.00
cor(index_actualactualhit,outcome_actual)                                          1.01
cor(index_counterfactualcounterfactualclose,outcome_actual)                        1.00
cor(index_counterfactualcounterfactualhit,outcome_actual)                          1.01
                                                                                   Bulk_ESS
sd(Intercept)                                                                          1772
sd(index_actualactualclose)                                                             910
sd(index_actualactualhit)                                                              1024
sd(index_counterfactualcounterfactualclose)                                            1102
sd(index_counterfactualcounterfactualhit)                                              1360
sd(outcome_actual)                                                                     1221
cor(Intercept,index_actualactualclose)                                                 3259
cor(Intercept,index_actualactualhit)                                                   2183
cor(index_actualactualclose,index_actualactualhit)                                     1812
cor(Intercept,index_counterfactualcounterfactualclose)                                 1238
cor(index_actualactualclose,index_counterfactualcounterfactualclose)                    311
cor(index_actualactualhit,index_counterfactualcounterfactualclose)                      507
cor(Intercept,index_counterfactualcounterfactualhit)                                   1487
cor(index_actualactualclose,index_counterfactualcounterfactualhit)                      330
cor(index_actualactualhit,index_counterfactualcounterfactualhit)                        363
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit)     1388
cor(Intercept,outcome_actual)                                                           881
cor(index_actualactualclose,outcome_actual)                                             556
cor(index_actualactualhit,outcome_actual)                                               689
cor(index_counterfactualcounterfactualclose,outcome_actual)                            1250
cor(index_counterfactualcounterfactualhit,outcome_actual)                              1664
                                                                                   Tail_ESS
sd(Intercept)                                                                          2323
sd(index_actualactualclose)                                                            1374
sd(index_actualactualhit)                                                              1620
sd(index_counterfactualcounterfactualclose)                                            1926
sd(index_counterfactualcounterfactualhit)                                              2216
sd(outcome_actual)                                                                     2105
cor(Intercept,index_actualactualclose)                                                 2034
cor(Intercept,index_actualactualhit)                                                   2671
cor(index_actualactualclose,index_actualactualhit)                                     2987
cor(Intercept,index_counterfactualcounterfactualclose)                                 2292
cor(index_actualactualclose,index_counterfactualcounterfactualclose)                    590
cor(index_actualactualhit,index_counterfactualcounterfactualclose)                     1166
cor(Intercept,index_counterfactualcounterfactualhit)                                   2044
cor(index_actualactualclose,index_counterfactualcounterfactualhit)                      784
cor(index_actualactualhit,index_counterfactualcounterfactualhit)                        872
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit)     2445
cor(Intercept,outcome_actual)                                                          1398
cor(index_actualactualclose,outcome_actual)                                            1563
cor(index_actualactualhit,outcome_actual)                                              1197
cor(index_counterfactualcounterfactualclose,outcome_actual)                            2261
cor(index_counterfactualcounterfactualhit,outcome_actual)                              2723

Population-Level Effects: 
                                                                      Estimate
Intercept                                                               -24.07
conditioncounterfactual_first                                            11.47
index_actualactualclose                                                   8.17
index_actualactualhit                                                    13.42
index_counterfactualcounterfactualclose                                 -30.20
index_counterfactualcounterfactualhit                                   -50.29
outcome_actual                                                           98.53
conditioncounterfactual_first:index_actualactualclose                    -5.39
conditioncounterfactual_first:index_actualactualhit                     -16.97
conditioncounterfactual_first:index_counterfactualcounterfactualclose    -8.13
conditioncounterfactual_first:index_counterfactualcounterfactualhit     -25.51
conditioncounterfactual_first:outcome_actual                             10.72
                                                                      Est.Error
Intercept                                                                  3.50
conditioncounterfactual_first                                              4.91
index_actualactualclose                                                    3.44
index_actualactualhit                                                      4.90
index_counterfactualcounterfactualclose                                    3.72
index_counterfactualcounterfactualhit                                      4.60
outcome_actual                                                             5.40
conditioncounterfactual_first:index_actualactualclose                      4.94
conditioncounterfactual_first:index_actualactualhit                        7.02
conditioncounterfactual_first:index_counterfactualcounterfactualclose      5.17
conditioncounterfactual_first:index_counterfactualcounterfactualhit        6.44
conditioncounterfactual_first:outcome_actual                               7.60
                                                                      l-95% CI
Intercept                                                               -30.87
conditioncounterfactual_first                                             1.62
index_actualactualclose                                                   1.47
index_actualactualhit                                                     3.74
index_counterfactualcounterfactualclose                                 -37.61
index_counterfactualcounterfactualhit                                   -59.57
outcome_actual                                                           87.81
conditioncounterfactual_first:index_actualactualclose                   -15.00
conditioncounterfactual_first:index_actualactualhit                     -30.40
conditioncounterfactual_first:index_counterfactualcounterfactualclose   -18.18
conditioncounterfactual_first:index_counterfactualcounterfactualhit     -38.16
conditioncounterfactual_first:outcome_actual                             -4.16
                                                                      u-95% CI
Intercept                                                               -17.24
conditioncounterfactual_first                                            21.22
index_actualactualclose                                                  14.79
index_actualactualhit                                                    22.69
index_counterfactualcounterfactualclose                                 -23.19
index_counterfactualcounterfactualhit                                   -41.47
outcome_actual                                                          108.98
conditioncounterfactual_first:index_actualactualclose                     4.35
conditioncounterfactual_first:index_actualactualhit                      -3.33
conditioncounterfactual_first:index_counterfactualcounterfactualclose     2.17
conditioncounterfactual_first:index_counterfactualcounterfactualhit     -12.98
conditioncounterfactual_first:outcome_actual                             26.02
                                                                      Rhat
Intercept                                                             1.00
conditioncounterfactual_first                                         1.00
index_actualactualclose                                               1.00
index_actualactualhit                                                 1.00
index_counterfactualcounterfactualclose                               1.00
index_counterfactualcounterfactualhit                                 1.00
outcome_actual                                                        1.00
conditioncounterfactual_first:index_actualactualclose                 1.00
conditioncounterfactual_first:index_actualactualhit                   1.00
conditioncounterfactual_first:index_counterfactualcounterfactualclose 1.00
conditioncounterfactual_first:index_counterfactualcounterfactualhit   1.00
conditioncounterfactual_first:outcome_actual                          1.00
                                                                      Bulk_ESS
Intercept                                                                 1592
conditioncounterfactual_first                                             1214
index_actualactualclose                                                   2313
index_actualactualhit                                                     1993
index_counterfactualcounterfactualclose                                   2457
index_counterfactualcounterfactualhit                                     1597
outcome_actual                                                            1983
conditioncounterfactual_first:index_actualactualclose                     2251
conditioncounterfactual_first:index_actualactualhit                       1990
conditioncounterfactual_first:index_counterfactualcounterfactualclose     2373
conditioncounterfactual_first:index_counterfactualcounterfactualhit       1212
conditioncounterfactual_first:outcome_actual                              1712
                                                                      Tail_ESS
Intercept                                                                 2398
conditioncounterfactual_first                                             1554
index_actualactualclose                                                   2743
index_actualactualhit                                                     2073
index_counterfactualcounterfactualclose                                   2375
index_counterfactualcounterfactualhit                                     1927
outcome_actual                                                            2818
conditioncounterfactual_first:index_actualactualclose                     2511
conditioncounterfactual_first:index_actualactualhit                       2024
conditioncounterfactual_first:index_counterfactualcounterfactualclose     2225
conditioncounterfactual_first:index_counterfactualcounterfactualhit       2082
conditioncounterfactual_first:outcome_actual                              2341

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    29.90      0.65    28.65    31.18 1.00     2704     3029

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

3.5.2.3 Specific hypotheses

df.exp2.causal.posterior = df.exp2.combined.long %>% 
  filter(question == "causal",
         clip <= 18) %>% 
  group_by(clip, index_actual, index_counterfactual, outcome_actual, condition) %>% 
  summarize(value = mean(rating)) %>% 
  ungroup() %>% 
  add_fitted_draws(newdata = .,
                   model = fit.brm_exp2_causal,
                   re_formula = NA) %>% 
  ungroup() %>% 
  select(clip, condition, .value, .draw) %>% 
  spread(clip, .value)

# difference between clips 13 and 14 versus 8 
df.exp2.causal.posterior %>% 
    mutate(difference = (`13`  + `14`) / 2 - `8`) %>% 
    pull(difference) %>% 
    mean_hdi() %>% 
    mutate_if(is.numeric, ~round(., 2)) %>% 
    summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])"))

# difference between clips 5 and 6 versus 11 
df.exp2.causal.posterior %>% 
    mutate(difference = (`5`  + `6`) / 2 - `11`) %>% 
    pull(difference) %>% 
    mean_hdi() %>% 
    mutate_if(is.numeric, ~round(., 2)) %>% 
    summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])"))

               difference
1 (-0.55 [-12.09, 10.97])
             difference
1 (-5.48 [-13.8, 2.73])

3.5.3 Tables

3.5.3.1 Bayesian regression

fit.brm_exp2_causal %>% 
  tidy() %>% 
  filter(str_detect(term, "b_")) %>% 
  mutate(term = str_remove(term, "b_"),
         term = tolower(term)) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  rename(`lower 95% HDI` = lower,
         `upper 95% HDI` = upper) %>% 
  xtable() %>%
  print(include.rownames = F,
        booktabs = T)

% latex table generated in R 3.6.1 by xtable 1.8-4 package
% Tue Nov 12 15:05:29 2019
\begin{table}[ht]
\centering
\begin{tabular}{lrrrr}
  \toprule
term & estimate & std.error & lower 95\% HDI & upper 95\% HDI \\ 
  \midrule
intercept & -24.07 & 3.50 & -29.80 & -18.27 \\ 
  conditioncounterfactual\_first & 11.47 & 4.91 & 3.29 & 19.40 \\ 
  index\_actualactualclose & 8.17 & 3.44 & 2.39 & 13.62 \\ 
  index\_actualactualhit & 13.42 & 4.90 & 5.42 & 21.30 \\ 
  index\_counterfactualcounterfactualclose & -30.20 & 3.72 & -36.41 & -24.36 \\ 
  index\_counterfactualcounterfactualhit & -50.29 & 4.60 & -57.94 & -42.84 \\ 
  outcome\_actual & 98.53 & 5.40 & 89.53 & 107.57 \\ 
  conditioncounterfactual\_first:index\_actualactualclose & -5.39 & 4.94 & -13.38 & 2.80 \\ 
  conditioncounterfactual\_first:index\_actualactualhit & -16.97 & 7.02 & -28.61 & -5.75 \\ 
  conditioncounterfactual\_first:index\_counterfactualcounterfactualclose & -8.13 & 5.17 & -16.51 & 0.35 \\ 
  conditioncounterfactual\_first:index\_counterfactualcounterfactualhit & -25.51 & 6.44 & -36.53 & -15.18 \\ 
  conditioncounterfactual\_first:outcome\_actual & 10.72 & 7.60 & -1.63 & 23.34 \\ 
   \bottomrule
\end{tabular}
\end{table}

4 Experiment 3: 3 balls

4.1 Clips

include_graphics("../../figures/diagrams/exp3_diagrams.png")

4.2 Read in data

# clipinfo
df.exp3.clipinfo = tibble(
    clip = 1:32,
    outcome_both = c(0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1),
    outcome_a = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1),
    outcome_b = c(0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1),
    outcome_none = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
    clipindex = rep(1:2, 16))

# COUNTERFACTUAL JUDGMENTS 
df.exp3.counterfactual = read.csv("../../data/experiment3_counterfactual.csv", header = T, stringsAsFactors = F)

# demographics
df.exp3.counterfactual.demographics = df.exp3.counterfactual$extra %>%
  strsplit("\n") %>%
  unlist() %>%
  strsplit(",") %>%
  unlist() %>%
  matrix(ncol = 6, byrow = T) %>%
  as.data.frame(stringsAsFactors = F) %>%
  mutate_all("as.numeric") %>%
  setNames(c("gender", "age", "difficulty", "smoothness", "time", "condition")) %>%
  mutate(
    participant = 1:nrow(.),
    condition = factor(condition, levels = 20:21, labels = c("A", "B")),
    gender = factor(gender, levels = 1:2, labels = c("female", "male"))
  ) %>% 
  select(participant, condition, gender, age, time, smoothness, difficulty)

# judgments
df.exp3.counterfactual.long = df.exp3.counterfactual$data %>%
  strsplit("\n") %>%
  unlist() %>%
  strsplit(",") %>%
  unlist() %>%
  matrix(ncol = 64, byrow = T) %>%
  as.data.frame(stringsAsFactors = F) %>%
  mutate_all("as.numeric") %>%
  setNames(str_c(c("clip_", "rating_"), rep(1:32, each = 2))) %>%
  mutate(participant = 1:nrow(.)) %>%
  gather("index", "value", -participant) %>% 
  separate(index, into = c("index", "order")) %>% 
  mutate(order = as.numeric(order)) %>% 
  spread(index, value) %>%
  left_join(select(df.exp3.counterfactual.demographics, participant, ball = condition),
            by = "participant") %>% 
  select(participant, ball, clip, everything()) %>%
  arrange(participant, clip)

# CAUSAL JUDGMENTS
df.exp3.causal = read.csv("../../data/experiment3_causal.csv", header = T, stringsAsFactors = F)

# demographics 
df.exp3.causal.demographics = df.exp3.causal$extra %>%
  strsplit("\n") %>%
  unlist() %>%
  strsplit(",") %>%
  unlist() %>%
  matrix(nrow = nrow(df.exp3.causal), byrow = T) %>%
  as.data.frame(stringsAsFactors = F) %>%
  mutate_all("as.numeric") %>%
  setNames(c("gender", "age", "difficulty", "smoothness", "time", "counterbalance", "replayType", "experiment")) %>%
  mutate(participant = 1:nrow(.),
         counterbalance = ifelse(participant %in% c(1, 3, 5, 6, 9), 1, counterbalance),
         gender = factor(gender, levels = 1:2, labels = c("female", "male"))) %>%
  select(participant, gender, age, counterbalance, time, smoothness, difficulty)
  
# judgments
df.exp3.causal.long = df.exp3.causal$data %>%
  strsplit("\n") %>%
  unlist() %>%
  strsplit(",") %>%
  unlist() %>%
  matrix(ncol = 192, byrow = T) %>%
  as.data.frame(stringsAsFactors = F) %>%
  mutate_all("as.numeric") %>%
  setNames(str_c(c("clip_", "X_", "Y_", "Z_", "rating1_", "rating2_"), rep(1:32, each = 6))) %>%
  select(contains("clip"), contains("rating")) %>%
  mutate(participant = 1:nrow(.)) %>%
  gather("index", "value", -participant) %>% 
  separate(index, into = c("index", "order")) %>% 
  mutate(order = as.numeric(order)) %>% 
  spread(index, value) %>% 
  left_join(select(df.exp3.causal.demographics, participant, counterbalance),
            by = "participant") %>% 
  arrange(participant, clip) %>% 
  mutate(
    A = ifelse(counterbalance == 1, rating1, rating2),
    B = ifelse(counterbalance == 1, rating2, rating1)
  ) %>% 
  select(-rating1, -rating2) %>% 
  gather("ball", "rating", c(A, B)) %>% 
  select(participant, clip, ball, rating, order) %>% 
  arrange(participant, clip, ball)

4.2.1 Demographics

# counterfactual condition 
df.exp3.counterfactual.demographics %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2))

# causal condition 
df.exp3.causal.demographics %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2))

   n age_mean age_sd n_female time_mean time_sd
1 80       33   10.1       34     18.08    4.63
   n age_mean age_sd n_female time_mean time_sd
1 41       34   10.5       21     21.19    4.96

4.3 Read in model predictions

# read model predictions 
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "3ball")]

df.exp3.model = files %>%
  map(~ read.csv(str_c(path, .))) %>% 
  reduce(rbind) %>%
  rename(clip = trial)

4.3.1 Counterfactual condition

# calculate mean counterfactual judgments 
df.exp3.counterfactual.means = df.exp3.counterfactual.long %>% 
  group_by(clip, ball) %>%
  summarize(rating_mean = mean(rating),
            rating_low = smean.cl.boot(rating)[2],
            rating_high = smean.cl.boot(rating)[3]) %>% 
  ungroup()

# find noisy simulation model that correlates best with the mean counterfactual judgments
df.exp3.model.counterfactual = df.exp3.model %>% 
  select(clip, contains("whether"), noise) %>% 
  gather("ball", "prediction", c(A_whether, B_whether)) %>% 
  mutate(ball = str_remove(ball, "_whether")) %>% 
  arrange(noise, clip, ball) %>% 
  left_join(df.exp3.clipinfo, by = "clip") %>% 
  mutate(prediction = ifelse(outcome_both == 1, 1 - prediction, prediction)) %>% 
  group_by(noise) %>% 
  nest() %>% 
  mutate(correlation = map(data, ~ cor(.x$prediction,
                                       df.exp3.counterfactual.means$rating_mean)),
         correlation = unlist(correlation)) %>% 
  ungroup() %>% 
  filter(correlation == max(correlation)) %>% 
  unnest(data) %>% 
  mutate(prediction = prediction * 100)

4.3.2 Causal condition

# calculate mean causal judgments 
df.exp3.causal.means = df.exp3.causal.long %>% 
  group_by(clip, ball) %>%
  summarize(rating_mean = mean(rating),
            rating_low = smean.cl.boot(rating)[2],
            rating_high = smean.cl.boot(rating)[3]) %>% 
  ungroup()

df.exp3.model.causal = df.exp3.model %>% 
  # take into account difference-making 
  mutate_at(.vars = vars(A_how:A_robust),
            .funs = list(~ . * A_difference)) %>% 
  mutate_at(.vars = vars(B_how:B_robust),
            .funs = list(~ . * B_difference)) %>% 
  # choose model based on best fit with counterfactual data 
  filter(noise == unique(df.exp3.model.counterfactual$noise)) %>%
  gather("index", "value", A_difference:B_robust) %>% 
  separate(index, into = c("ball", "aspect")) %>% 
  spread(aspect, value) %>% 
  select(clip, ball:whether)

4.3.3 Heuristic model

• clean up this section (use map instead of for loop)

l.features = fromJSON("data/features.json")

df.features.initial = l.features[["appearance"]] %>%
  cbind(l.features[["initialVelocity"]][["ballA"]]) %>%
  cbind(l.features[["initialVelocity"]][["ballB"]]) %>%
  cbind(l.features[["initialVelocity"]][["ballE"]]) %>%
  setNames(c(
    paste0("ball", LETTERS[c(1, 2, 5)], "_appearance"),
    "ballA_velx",
    "ballA_vely",
    "ballB_velx",
    "ballB_vely",
    "ballE_velx",
    "ballE_vely"
  )) %>%
  mutate(clip = 1:nrow(.)) %>%
  pivot_longer(cols = starts_with("ball")) %>%
  separate(name, into = c("ball", "property")) %>%
  pivot_wider(names_from = property, values_from = value) %>% 
  mutate(ball = str_remove_all(ball, pattern = "ball"))

# information about collisions
df.features.collisions = data.frame()
ballnames = c("ballA", "ballB", "ballE")

for (i in 1:32) {
  df.tmp = data.frame()
  for (j in 1:length(ballnames)) {
    ncollisions = length(l.features[["collisions"]][[ballnames[j]]][["object"]][[i]])
    if (ncollisions > 0) {
      tmp = data.frame(
        ball = ballnames[j],
        object = l.features[["collisions"]][[ballnames[j]]][["object"]][[i]],
        time = l.features[["collisions"]][[ballnames[j]]][["time"]][[i]],
        pre.velx = l.features[["collisions"]][[ballnames[j]]][["pre"]][[i]][["x"]],
        pre.vely = l.features[["collisions"]][[ballnames[j]]][["pre"]][[i]][["y"]],
        post.velx = l.features[["collisions"]][[ballnames[j]]][["post"]][[i]][["x"]],
        post.vely = l.features[["collisions"]][[ballnames[j]]][["post"]][[i]][["y"]],
        ncollision = 1:ncollisions
      )
      df.tmp = df.tmp %>%
        rbind(tmp)
    }
  }
  df.features.collisions = df.features.collisions %>%
    rbind(df.tmp %>%
      mutate(clip = i) %>%
      select(clip, ball, everything()))
}

# find end of clip
tmp = df.features.collisions %>%
  group_by(clip) %>%
  filter(
    ball == "ballE",
    str_detect(object, "goal|Left|Right")
  ) %>%
  group_by(clip) %>%
  mutate(endclip = max(time)) %>%
  select(clip, endclip) %>%
  ungroup() %>%
  distinct()

# remove events that happen after the end of the clip
df.features.collisions = df.features.collisions %>%
  left_join(tmp) %>%
  mutate(endclip = ifelse(is.na(endclip), 400, endclip)) %>%
  group_by(clip) %>%
  filter(time <= endclip)

Joining, by = "clip"

# E initially at rest or moving?
tmp.movingE = df.features.initial %>%
  filter(ball == "E") %>%
  mutate(E.moving = ifelse(velx == 0 & vely == 0, 1, 0)) %>%
  select(clip, E.moving)

# contact with ball E
tmp.contactE = df.features.collisions %>%
  filter(object == "ballE") %>%
  mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
  mutate(contactE = 1) %>%
  select(clip, ball, contactE) %>%
  group_by(clip, ball) %>%
  summarise(contactE = any(contactE %>% as.logical()) * 1) %>%
  left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .) %>%
  mutate(contactE = ifelse(is.na(contactE), 0, contactE))

Joining, by = c("clip", "ball")

# number of collisions
tmp.ncollisions = df.features.collisions %>%
  filter(
    str_detect(object, "ball"),
    ball != "ballE"
  ) %>%
  mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
  group_by(clip, ball) %>%
  count() %>%
  rename(ncollision = n)

# exclusivity (no contact with ball E by the other ball)
tmp.exclusive = df.features.collisions %>%
  filter(
    str_detect(object, "ball"),
    ball != "ballE"
  ) %>%
  count(clip, ball, object) %>%
  mutate(
    AE = ifelse(ball == "ballA" & object == "ballE", 1, 0),
    BE = ifelse(ball == "ballB" & object == "ballE", 1, 0)
  ) %>%
  group_by(clip) %>%
  summarise(
    AE = any(AE %>% as.logical()) * 1,
    BE = any(BE %>% as.logical()) * 1
  ) %>%
  mutate(
    A = ifelse(AE == 1 & BE == 0, 1, 0),
    B = ifelse(AE == 0 & BE == 1, 1, 0)
  ) %>%
  select(clip, A, B) %>%
  gather(ball, exclusive, -clip)

# impact
func_angle = function(x, y) {
  dot.prod = x %*% y
  norm.x = norm(x, type = "2")
  norm.y = norm(y, type = "2")
  theta = acos(dot.prod / (norm.x * norm.y))
  as.numeric(theta)
}

# impact on ball E 
tmp.impactE = df.features.collisions %>%
  rowwise() %>%
  filter(
    str_detect(object, "ball"),
    ball == "ballE"
  ) %>%
  mutate(
    pre.speed = abs(pre.velx) + abs(pre.vely),
    post.speed = abs(post.velx) + abs(post.vely),
    speed.diff = post.speed - pre.speed,
    direction.diff = func_angle(c(pre.velx, pre.vely), c(post.velx, post.vely)),
    direction.diff = ifelse(is.na(direction.diff), abs(atan2(post.vely - pre.vely, post.velx - pre.velx)), direction.diff)
  ) %>%
  ungroup() %>%
  group_by(clip, object) %>%
  summarise(
    speed.diff = sum(speed.diff),
    direction.diff = sum(direction.diff)
  ) %>%
  rename(ball = object) %>%
  mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
  left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .) %>%
  mutate_at(.vars = vars(speed.diff, direction.diff), .funs = ~ ifelse(is.na(.), 0, .)) %>%
  arrange(clip, ball) %>%
  rename(
    E.speed.diff = speed.diff,
    E.direction.diff = direction.diff
  )

Joining, by = c("clip", "ball")

# total impact 
tmp.impactTotal = df.features.collisions %>%
  rowwise() %>%
  filter(str_detect(object, "ballA|ballB")) %>%
  mutate(
    pre.speed = abs(pre.velx) + abs(pre.vely),
    post.speed = abs(post.velx) + abs(post.vely),
    speed.diff = post.speed - pre.speed,
    direction.diff = func_angle(c(pre.velx, pre.vely), c(post.velx, post.vely)),
    direction.diff = ifelse(is.na(direction.diff), abs(atan2(post.vely - pre.vely, post.velx - pre.velx)), direction.diff)
  ) %>%
  ungroup() %>%
  group_by(clip, object) %>%
  summarise(
    speed.diff = sum(speed.diff),
    direction.diff = sum(direction.diff)
  ) %>%
  rename(ball = object) %>%
  mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
  left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .) %>%
  mutate_at(.vars = vars(speed.diff, direction.diff), .funs = ~ ifelse(is.na(.), 0, .)) %>%
  arrange(clip, ball) %>%
  rename(
    total.speed.diff = speed.diff,
    total.direction.diff = direction.diff
  )

Joining, by = c("clip", "ball")

# transfer of force 
tmp.transfer = df.features.collisions %>%
  filter(
    str_detect(ball, "ballA|ballB"),
    str_detect(object, "ball")
  ) %>%
  mutate(
    AE = ifelse(ball == "ballA" & object == "ballE", 1, NA),
    BE = ifelse(ball == "ballB" & object == "ballE", 1, NA),
    AB = ifelse((ball == "ballA" & object == "ballB"), 1, NA)
  ) %>%
  mutate_at(.vars = vars(AE, BE, AB), .funs = ~ . * time) %>%
  filter(!is.na(AE) | !is.na(BE) | !is.na(AB)) %>%
  arrange(clip, time) %>%
  group_by(clip) %>%
  summarise(
    A = any(!is.na(AE)),
    B = any(!is.na(BE)),
    A = ifelse(any(!is.na(AB)) & any(!is.na(BE)) & max(AB, na.rm = T) < min(BE, na.rm = T), T, A),
    B = ifelse(any(!is.na(AB)) & any(!is.na(AE)) & max(AB, na.rm = T) < min(AE, na.rm = T), T, B)
  ) %>%
  gather(ball, transfer, -clip) %>%
  arrange(clip, ball)

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf

Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf

Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf

Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf

Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in min(AE, na.rm = T): no non-missing arguments to min; returning
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in min(AE, na.rm = T): no non-missing arguments to min; returning
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in min(AE, na.rm = T): no non-missing arguments to min; returning
Inf

Warning in min(AE, na.rm = T): no non-missing arguments to min; returning
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in min(AE, na.rm = T): no non-missing arguments to min; returning
Inf

Warning in min(AE, na.rm = T): no non-missing arguments to min; returning
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in min(AE, na.rm = T): no non-missing arguments to min; returning
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in min(AE, na.rm = T): no non-missing arguments to min; returning
Inf

Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf

Warning in min(AE, na.rm = T): no non-missing arguments to min; returning
Inf

# collect features
df.features = df.features.initial %>%
  filter(ball != "E") %>%
  mutate(ball = factor(ball)) %>%
  mutate(
    moving = ifelse(velx == 0 & vely == 0, 0, 1),
    speed = abs(velx) + abs(vely)
  ) %>%
  left_join(tmp.contactE) %>%
  left_join(tmp.impactE) %>%
  left_join(tmp.impactTotal) %>%
  left_join(tmp.transfer %>% mutate(transfer = transfer * 1)) %>%
  left_join(tmp.movingE) %>%
  left_join(tmp.exclusive) %>%
  select(-c(appearance, velx, vely))

Joining, by = c("clip", "ball")

Joining, by = c("clip", "ball")
Joining, by = c("clip", "ball")
Joining, by = c("clip", "ball")

Warning: Column `ball` joining factor and character vector, coercing into
character vector

Joining, by = "clip"
Joining, by = c("clip", "ball")

# remove temporary variables
rm(list = ls()[str_detect(ls(), "tmp")])

4.4 Counterfactual condition

4.4.1 Plots

4.4.1.1 Bar graph

df.plot = df.exp3.counterfactual.long %>%
  left_join(df.exp3.model.counterfactual %>%
    select(clip, ball, prediction) %>%
    mutate(ball = as.factor(ball)),
  by = c("clip", "ball")
  ) %>%
  left_join(df.exp3.clipinfo,
            by = "clip") %>% 
  gather("model", "value", c(rating, prediction)) %>%
  mutate(
    colorindex = 1,
    colorindex = ifelse(model == "rating" & ball == "A" & outcome_b == 1, 2, colorindex),
    colorindex = ifelse(model == "rating" & ball == "B" & outcome_a == 1, 2, colorindex),
    colorindex = ifelse(model == "prediction", 3, colorindex),
    colorindex = as.factor(colorindex),
    model = factor(model, levels = c("rating", "prediction"))
  ) %>%
  arrange(participant, clip, ball)

df.text = df.plot %>%
  expand(ball, clip, model) %>%
  mutate(
    label = ifelse(model == "rating" & ball == "A", clip, NA),
    y = 110,
    x = 1.7,
    colorindex = NA
  )

actual_counterfactual_threeballs_plot(df.plot, "counterfactual judgment")
# ggsave("../../figures/plots/exp3_counterfactual_bars.pdf",
#        width = 12,
#        height = 6)

4.4.2 Stats

4.4.2.1 Model evaluation

# best fitting model 
df.exp3.counterfactual.means %>% 
  left_join(df.exp3.model.counterfactual) %>% 
  summarize(simulation_r = cor(rating_mean, prediction),
            simulation_rmse = rmse(rating_mean, prediction)) %>% 
  mutate_all(~ round(., 2))

Joining, by = c("clip", "ball")

Warning: Column `ball` joining factor and character vector, coercing into
character vector

# deterministic model 
df.exp3.counterfactual.means %>% 
  left_join(df.exp3.clipinfo %>% 
              select(clip, outcome_a, outcome_b) %>% 
              gather("ball", "outcome", -clip) %>% 
              mutate(ball = factor(ball,
                                   levels = c("outcome_a", "outcome_b"),
                                   labels = c("B", "A")))) %>% 
  mutate(outcome = outcome * 100) %>% 
  summarize(simulation_r = cor(rating_mean, outcome),
            simulation_rmse = rmse(rating_mean, outcome)) %>% 
  mutate_all(~ round(., 2))

Joining, by = c("clip", "ball")

Warning: Column `ball` joining factors with different levels, coercing to
character vector

# A tibble: 1 x 2
  simulation_r simulation_rmse
         <dbl>           <dbl>
1         0.87            20.7
# A tibble: 1 x 2
  simulation_r simulation_rmse
         <dbl>           <dbl>
1         0.83            30.3

4.5 Causal condition

4.5.1 Stats

4.5.1.1 Bayesian mixed effects model

df.data = df.exp3.causal.long %>% 
  left_join(df.exp3.model.causal,
            by = c("clip", "ball"))

fit.brm_exp3_causal_w = brm(formula = rating ~ 1 + whether + (1 + whether | participant),
                                  data = df.data,
                                  save_all_pars = T,
                                  cores = 2,
                                  seed = 1,
                                  file = "cache/brm_exp3_causal_w")

fit.brm_exp3_causal_wh = brm(formula = rating ~ 1 + whether + how + (1 + whether + how | participant),
                                   data = df.data,
                                   save_all_pars = T,
                                   cores = 2,
                                   seed = 1,
                                   file = "cache/brm_exp3_causal_wh")

fit.brm_exp3_causal_whs = brm(formula = rating ~ 1 + whether + how + sufficient + (1 + whether + how + sufficient | participant),
                                    data = df.data,
                                    save_all_pars = T,
                                    cores = 2,
                                    seed = 1,
                                    file = "cache/brm_exp3_causal_whs")

fit.brm_exp3_causal_w %>% summary()
fit.brm_exp3_causal_wh %>% summary()
fit.brm_exp3_causal_whs %>% summary()

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: rating ~ 1 + whether + (1 + whether | participant) 
   Data: df.data (Number of observations: 2624) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~participant (Number of levels: 41) 
                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(Intercept)              7.76      1.64     4.78    11.03 1.00     1762
sd(whether)                8.39      2.93     2.47    13.90 1.01      677
cor(Intercept,whether)     0.00      0.37    -0.59     0.81 1.00     1100
                       Tail_ESS
sd(Intercept)              2343
sd(whether)                 995
cor(Intercept,whether)     1344

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    31.68      1.61    28.41    34.87 1.00     2557     2982
whether      39.30      2.16    35.11    43.54 1.00     3175     2876

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    35.64      0.51    34.65    36.66 1.00     4587     2904

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: rating ~ 1 + whether + how + (1 + whether + how | participant) 
   Data: df.data (Number of observations: 2624) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~participant (Number of levels: 41) 
                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(Intercept)              5.27      2.01     1.24     9.36 1.00     1414
sd(whether)               13.10      2.41     8.76    18.05 1.00     1931
sd(how)                   13.09      2.19     9.24    17.93 1.00     1141
cor(Intercept,whether)     0.29      0.32    -0.34     0.86 1.00      624
cor(Intercept,how)        -0.18      0.32    -0.69     0.53 1.00      467
cor(whether,how)          -0.63      0.16    -0.87    -0.25 1.00      757
                       Tail_ESS
sd(Intercept)              1535
sd(whether)                2994
sd(how)                    1965
cor(Intercept,whether)      987
cor(Intercept,how)          640
cor(whether,how)           1864

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    10.29      1.56     7.30    13.41 1.00     3821     3177
whether      28.40      2.68    23.18    33.56 1.00     2444     3161
how          36.07      2.60    31.01    41.17 1.00     2066     2624

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    31.58      0.44    30.72    32.48 1.00     6371     3175

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: rating ~ 1 + whether + how + sufficient + (1 + whether + how + sufficient | participant) 
   Data: df.data (Number of observations: 2624) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~participant (Number of levels: 41) 
                          Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                 6.27      2.08     1.93    10.42 1.00
sd(whether)                  13.57      2.22     9.65    18.37 1.00
sd(how)                      16.14      2.43    11.89    21.26 1.01
sd(sufficient)               14.49      2.84     9.35    20.52 1.00
cor(Intercept,whether)        0.12      0.27    -0.38     0.66 1.01
cor(Intercept,how)           -0.16      0.26    -0.62     0.38 1.01
cor(whether,how)             -0.71      0.13    -0.90    -0.42 1.01
cor(Intercept,sufficient)    -0.01      0.28    -0.54     0.58 1.00
cor(whether,sufficient)       0.86      0.10     0.62     0.98 1.00
cor(how,sufficient)          -0.68      0.15    -0.91    -0.34 1.00
                          Bulk_ESS Tail_ESS
sd(Intercept)                  868      589
sd(whether)                   1723     2357
sd(how)                       1117     2189
sd(sufficient)                2004     2611
cor(Intercept,whether)         464      668
cor(Intercept,how)             438      794
cor(whether,how)               495     1618
cor(Intercept,sufficient)      565      845
cor(whether,sufficient)       1779     2769
cor(how,sufficient)           1543     2706

Population-Level Effects: 
           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     10.43      1.59     7.30    13.50 1.00     3044     3001
whether       27.93      2.60    22.80    33.09 1.00     1382     2074
how           26.47      2.91    20.83    32.42 1.00     1248     1933
sufficient    30.90      3.06    25.03    36.93 1.00     1643     2977

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    30.00      0.43    29.18    30.86 1.00     4986     3192

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

4.5.1.2 Model comparison

fit.brm_exp3_causal_w = add_criterion(fit.brm_exp3_causal_w,
                                            criterion = c("loo", "waic"),
                                            reloo = T)

No problematic observations found. Returning the original 'loo' object.

fit.brm_exp3_causal_wh = add_criterion(fit.brm_exp3_causal_wh,
                                             criterion = c("loo", "waic"),
                                             reloo = T)

No problematic observations found. Returning the original 'loo' object.

fit.brm_exp3_causal_whs = add_criterion(fit.brm_exp3_causal_whs,
                                              criterion = c("loo", "waic"),
                                              reloo = T)

No problematic observations found. Returning the original 'loo' object.

loo_compare(fit.brm_exp3_causal_w,
            fit.brm_exp3_causal_wh,
            fit.brm_exp3_causal_whs)

# loo_compare(fit.brm_exp3_causal_w,
#             fit.brm_exp3_causal_wh)

                        elpd_diff se_diff
fit.brm_exp3_causal_whs    0.0       0.0 
fit.brm_exp3_causal_wh  -125.2      15.6 
fit.brm_exp3_causal_w   -426.2      31.1 

4.5.1.3 Heuristic model

• z-scored predictors and invidiual responses

# Regression based on features
df.regression.features = df.exp3.causal.long %>%
  left_join(df.exp3.clipinfo, by = "clip") %>% 
  left_join(df.features %>% 
              mutate_at(.vars = vars(moving:exclusive), .funs = ~ scale(.)[,]),
            by = c("clip", "ball")) %>%
  clean_names() %>% 
  mutate(e_moving = 1 - e_moving)

# restrict the weights such that all predictors are positive
priors = c(set_prior("normal(0,10)", class = "b", lb = 0))

fit.brm_exp3_causal_heuristic = brm(formula = rating ~ moving + 
                                      speed + 
                                      contact_e + 
                                      e_speed_diff + 
                                      e_direction_diff + 
                                      total_speed_diff + 
                                      total_direction_diff + 
                                      transfer + 
                                      e_moving + 
                                      exclusive + (1 | participant),
                                    prior = priors,
                                    data = df.regression.features %>% 
                                      select(-c(clip, ball, order, clipindex, contains("outcome"))),
                                    save_all_pars = T,
                                    cores = 2,
                                    seed = 1,
                                    file = "cache/brm_exp3_causal_heuristic")
fit.brm_exp3_causal_heuristic

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: rating ~ moving + speed + contact_e + e_speed_diff + e_direction_diff + total_speed_diff + total_direction_diff + transfer + e_moving + exclusive + (1 | participant) 
   Data: df.regression.features %>% select(-c(clip, ball, o (Number of observations: 2624) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~participant (Number of levels: 41) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     8.53      1.22     6.44    11.18 1.00     1163     1936

Population-Level Effects: 
                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept               49.79      1.46    46.91    52.74 1.00     1129
moving                   0.22      0.22     0.01     0.81 1.00     3399
speed                    2.06      0.85     0.44     3.67 1.00     2319
contact_e                0.38      0.35     0.01     1.29 1.00     3201
e_speed_diff             0.12      0.12     0.00     0.46 1.00     4435
e_direction_diff         1.04      0.74     0.04     2.71 1.00     2202
total_speed_diff         2.21      0.95     0.41     4.17 1.00     2208
total_direction_diff     3.99      0.91     2.14     5.68 1.00     2718
transfer                15.59      0.79    14.06    17.14 1.00     3555
e_moving                 0.18      0.17     0.01     0.62 1.00     3296
exclusive                4.39      0.73     2.97     5.81 1.00     3526
                     Tail_ESS
Intercept                1854
moving                   1547
speed                    1266
contact_e                2020
e_speed_diff             2268
e_direction_diff         1685
total_speed_diff         1393
total_direction_diff     1742
transfer                 2992
e_moving                 1738
exclusive                2538

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    33.21      0.46    32.33    34.10 1.00     5313     2845

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

4.5.1.4 Predictions for mean causal judgments

df.exp3.causal.regression = df.exp3.causal.means %>% 
  mutate(w = fit.brm_exp3_causal_w %>% 
           fitted(newdata = df.exp3.model.causal,
                  re_formula = NA) %>% 
           as_tibble() %>% 
           pull(Estimate),
         wh = fit.brm_exp3_causal_wh %>% 
           fitted(newdata = df.exp3.model.causal,
                  re_formula = NA) %>% 
           as_tibble() %>% 
           pull(Estimate),
         whs = fit.brm_exp3_causal_whs %>% 
           fitted(newdata = df.exp3.model.causal,
                  re_formula = NA) %>% 
           as_tibble() %>% 
           pull(Estimate),
         heuristic = fit.brm_exp3_causal_heuristic %>% 
           fitted(newdata = df.exp3.model.causal %>% 
                    left_join(df.features, by = c("clip", "ball")) %>% 
                    clean_names(),
                  re_formula = NA) %>% 
           as_tibble() %>% 
           pull(Estimate))

4.5.1.5 Individual participant regressions

df.fit = df.exp3.causal.long %>% 
  left_join(df.exp3.model.causal,
            by = c("clip", "ball"))

# priors 
priors = c(set_prior("normal(0,10)", class = "b", lb = 0))

# initial model fits (used for compilation)
fit.brm_exp3_causal_individual_baseline = brm(formula = rating ~ 1,
                                              data = df.fit %>% 
                                                filter(participant == 1),
                                              save_all_pars = T,
                                              cores = 2,
                                              seed = 1,
                                              file = str_c("cache/brm_exp3_causal_individual_baseline"))

fit.brm_exp3_causal_individual_w = brm(formula = rating ~ 1 + whether,
                                       data = df.fit %>% 
                                         filter(participant == 1),
                                       prior = priors,
                                       save_all_pars = T,
                                       cores = 2,
                                       seed = 1,
                                       file = str_c("cache/brm_exp3_causal_individual_w"))

fit.brm_exp3_causal_individual_wh = brm(formula = rating ~ 1 + whether + how,
                                        data = df.fit %>% 
                                          filter(participant == 1),
                                        prior = priors,
                                        save_all_pars = T,
                                        cores = 2,
                                        seed = 1,
                                        file = str_c("cache/brm_exp3_causal_individual_wh"))

fit.brm_exp3_causal_individual_whs = brm(formula = rating ~ 1 + whether + how + sufficient,
                                         data = df.fit %>% 
                                           filter(participant == 1),
                                         prior = priors,
                                         save_all_pars = T,
                                         cores = 2,
                                         seed = 1,
                                         file = str_c("cache/brm_exp3_causal_individual_whs"))

# update model fits for different participants
df.exp3.causal.individual_fit = df.fit %>%
  group_by(participant) %>%
  nest() %>%
  ungroup() %>% 
  # fit model for each participant
  mutate(fit_baseline = map(.x = data,
                            .f = ~ update(fit.brm_exp3_causal_individual_baseline,
                                          newdata = .x)),
         fit_w = map(.x = data,
                     .f = ~ update(fit.brm_exp3_causal_individual_w,
                                   newdata = .x)),
         fit_wh = map(.x = data,
                      .f = ~ update(fit.brm_exp3_causal_individual_wh,
                                    newdata = .x)),
         fit_whs = map(.x = data,
                       .f = ~ update(fit.brm_exp3_causal_individual_whs,
                                     newdata = .x))) %>% 
  mutate(fit_baseline = map(.x = fit_baseline,
                            ~ add_criterion(.x, criterion = c("loo", "waic"))),
         fit_w = map(.x = fit_w,
                     ~ add_criterion(.x, criterion = c("loo", "waic"))),
         fit_wh = map(.x = fit_wh,
                      ~ add_criterion(.x, criterion = c("loo", "waic"))),
         fit_whs = map(.x = fit_whs,
                       ~ add_criterion(.x, criterion = c("loo", "waic"))),
         r_w = map2_dbl(.x = data,
                        .y = fit_w,
                        .f = ~ cor(.x$rating, .y %>% 
                                     fitted() %>% 
                                     .[, 1])),
         r_wh = map2_dbl(.x = data,
                         .y = fit_wh,
                         .f = ~ cor(.x$rating, .y %>% 
                                      fitted() %>% 
                                      .[, 1])),
         r_whs = map2_dbl(.x = data,
                          .y = fit_whs,
                          .f = ~ cor(.x$rating, .y %>% 
                                       fitted() %>% 
                                       .[, 1])),
         rmse_w = map2_dbl(.x = data,
                        .y = fit_w,
                        .f = ~ rmse(.x$rating, .y %>% 
                                     fitted() %>% 
                                     .[, 1])),
         rmse_wh = map2_dbl(.x = data,
                         .y = fit_wh,
                         .f = ~ rmse(.x$rating, .y %>% 
                                      fitted() %>% 
                                      .[, 1])),
         rmse_whs = map2_dbl(.x = data,
                          .y = fit_whs,
                          .f = ~ rmse(.x$rating, .y %>% 
                                       fitted() %>% 
                                       .[, 1])),
         model_comparison = pmap(.l = list(baseline = fit_baseline,
                                           w = fit_w,
                                           wh = fit_wh,
                                           whs = fit_whs),
                                 .f = ~ loo_compare(..1, ..2, ..3, ..4)),
         best_model = map_chr(.x = model_comparison,
                              .f = ~ rownames(.) %>% .[1]),
         best_model = factor(best_model,
                             levels = c("..1", "..2", "..3", "..4"),
                             labels = c("baseline", "w", "wh", "whs")))

save(list = c("df.exp3.causal.individual_fit"),
     file = "cache/exp3_causal_individual_fit.RData")

4.5.1.6 Load individual participant regressions

load("cache/exp3_causal_individual_fit.RData")

# count how many participants are best fit by the different models
df.exp3.causal.individual_fit %>% 
  count(best_model)

# A tibble: 2 x 2
  best_model     n
  <fct>      <int>
1 wh             2
2 whs           39

4.5.2 Plots

4.5.2.1 Bar graph

model_index = "whs"

# model predictions 
model_prediction = list(fit.brm_exp3_causal_w,
                        fit.brm_exp3_causal_wh,
                        fit.brm_exp3_causal_whs) %>%
  map(~ fitted(., df.exp3.causal.means %>% 
                 left_join(df.exp3.model.causal,
                           by = c("clip", "ball")),
               re_formula = NA)) %>% 
  reduce(rbind) %>% 
  as_tibble() %>% 
  mutate(ball = rep(c("A", "B"), n()/2),
         clip = rep(rep(1:32, each = 2), 3),
         model = rep(c("w", "wh", "whs"), each = 64))

df.plot = df.exp3.causal.long %>%
  left_join(df.exp3.clipinfo %>% 
              select(clip, outcome_both),
            by = c("clip")) %>% 
  left_join(model_prediction %>% 
              filter(model == model_index) %>% 
              select(model = Estimate,
                     clip,
                     ball),
            by = c("clip", "ball")) %>% 
  gather("model","value",c(rating, model)) %>%
    mutate(
    colorindex = 1,
    colorindex = ifelse(model == "rating" & outcome_both == 0, 1, colorindex),
    colorindex = ifelse(model == "rating" & outcome_both == 1, 2, colorindex),
    colorindex = ifelse(model == "model", 3, colorindex),
    colorindex = as.factor(colorindex),
    model = factor(model, levels = c("rating", "model"))
  ) %>%
  arrange(participant, clip, ball)

df.text = df.plot %>%
  expand(ball, clip, model) %>% 
  mutate(label = ifelse(model == "rating" & ball == "A", clip, NA),
         y = 110,
         x = 1.7,
         colorindex = NA
  )

actual_counterfactual_threeballs_plot(df.plot,"causal responsibility")
# ggsave("../../figures/plots/exp3_causal_bars.pdf",
#        width = 12,
#        height = 6)

4.5.2.2 Selection

# \u2713 = check mark 
# \u2717 = cross mark 

check = "\u2713"
cross = "\u2717"

x_labels = c(str_c("\nW ", cross, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross))
  

df.plot = df.exp3.causal.long %>%
  filter(clip %in% c(3, 7, 15, 16, 23)) %>%
  group_by(clip, ball) %>%
  summarise(mean = mean(rating),
            low = smean.cl.boot(rating)[2],
            high = smean.cl.boot(rating)[3]) %>%
  left_join(df.exp3.causal.regression %>% select(clip, ball, w, wh, whs),
            by = c("clip", "ball")) %>%
  ungroup() %>%
  gather("index", "value", c(mean, w, wh, whs)) %>%
  mutate_at(.vars = vars(low, high), .funs = ~ ifelse(index == "mean", ., NA)) %>%
  mutate(index = factor(index, levels = c("mean", "w", "wh", "whs")),
         clip = factor(clip, levels = c(7, 23, 3, 15, 16),
                       labels = c("causal chain", "double prevention", "joint causation",
                                  "overdetermination", "preemption")))

df.labels = df.plot %>% 
  distinct(clip, ball) %>% 
  arrange(clip, ball) %>% 
  mutate(labels = x_labels,
         value = -10,
         index = NA)

func_load_image = function(clip){
  readPNG(str_c("../../figures/diagrams/exp3_clip", clip, ".png"))
}

df.clips = df.plot %>% 
  distinct(clip) %>% 
  arrange(clip) %>% 
  mutate(number = c(7, 23, 3, 15, 16),
         grob = map(.x = number, .f = ~ func_load_image(clip = .x)),
         index = NA,
         value = NA,
         ball = NA,
         label = str_c("clip ", number))

ball_a = readPNG("../../figures/diagrams/ballA.png")
ball_a = rasterGrob(ball_a, interpolate = TRUE)

ball_b = readPNG("../../figures/diagrams/ballB.png")
ball_b = rasterGrob(ball_b, interpolate = TRUE)

p = ggplot(data = df.plot, 
       mapping = aes(x = ball,
                     y = value,
                     group = index,
                     fill = index)) +
  geom_bar(stat = "identity", color = "black",
           position = position_dodge(0.9),
           width = 0.9) +
  geom_errorbar(mapping = aes(ymin = low, ymax = high),
                width = 0,
                size = 1,
                position = position_dodge(0.9)) +
  annotation_custom(grob = ball_a, xmin = 0.5, xmax = 1.5, ymin = -25, ymax = -7) + 
  annotation_custom(grob = ball_b, xmin = 1.5, xmax = 2.5, ymin = -25, ymax = -7) +
  geom_text(data = df.labels,
            mapping = aes(label = labels,
                          y = -Inf),
            vjust = 1.2,
            family = "Arial Unicode MS",
            size = 8) + 
  geom_custom(data = df.clips,
              mapping = aes(data = grob, x = 1.5, y = Inf),
              grob_fun = function(x) rasterGrob(x,
                                                interpolate = T,
                                                vjust = -0.25)) + 
  geom_text(data = df.clips,
            mapping = aes(label = label,
                          y = Inf,
                          x = -Inf),
            size = 9,
            hjust = -0.2,
            vjust = -4) + 
  facet_grid(~clip) +
  scale_fill_grey(start = 1,
                  end = 0,
                  labels = c("mean rating  ", expression(CSM[W] ~ " "),
                             expression(CSM[WH] ~ " "),
                             expression(CSM[WHS] ~ " "))) +
  labs(y = "causal responsibility", fill = "") +
  coord_cartesian(clip = "off") + 
  theme_bw() +
  theme(legend.text.align = 0,
        text = element_text(size = 20),
        panel.grid = element_blank(),
        legend.position = "bottom",
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        legend.text = element_text(size = 30),
        strip.text = element_text(size = 30),
        legend.spacing = unit(0.5, "cm"),
        legend.background = element_rect(fill = "transparent"),
        legend.margin = margin(t = 5, unit = "cm"),
        plot.margin = margin(t = 8, l = 0.2, r = 0.2, b = 0.1, unit = "cm"))

# ggsave("../../figures/plots/exp3_selection_bars.pdf",
#        width = 20,
#        height = 10,
#        device = cairo_pdf)

4.5.2.3 Scatter plot

model = "w"
model = "wh"
model = "whs"
model = "heuristic"


if(model == "heuristic"){
  xlabel = "Heuristic"
}else{
  xlabel = bquote(CSM[.(toupper(model))])
}

l.models = list(w = fit.brm_exp3_causal_w, 
                wh = fit.brm_exp3_causal_wh,
                whs = fit.brm_exp3_causal_whs,
                heuristic = fit.brm_exp3_causal_heuristic)

df.data = df.exp3.causal.means %>% 
  left_join(df.exp3.model.causal, by = c("clip", "ball")) %>% 
  left_join(df.regression.features, by = c("clip", "ball"))

df.plot = l.models[[model]] %>%
  fitted(newdata = df.data,
         re_formula = NA) %>% 
  as_tibble() %>% 
  clean_names() %>%
  bind_cols(df.data)

ggplot(data = df.plot,
       mapping = aes(x = estimate, 
                     y = rating_mean)) +
  geom_abline(intercept = 0,
              slope = 1,
              linetype = 2) + 
  geom_smooth(aes(y = estimate,
                  ymin = q2_5,
                  ymax = q97_5),
              stat = "identity",
              color = "black") +
  geom_linerange(size = 0.75, 
                 mapping = aes(ymin = rating_low,
                               ymax = rating_high),
                 color = "gray50") +
  geom_point(size = 2) + 
  scale_color_manual(values = c("black", "#e41a1c"),
                     guide = F) +
  coord_fixed(xlim = c(0, 100),
              ylim = c(0, 100)) +
  labs(x = xlabel,
       y = "responsibility judgment") +
  annotate(geom = "text",
           label = str_c(
             "r = ", cor(df.plot$estimate, df.plot$rating_mean) %>% round(2) %>% as.character() %>% str_sub(start = 2),
             "\nRMSE = ", rmse(df.plot$estimate, df.plot$rating_mean) %>% round(2)),
           x = 5,
           y = 95,
           size = 6,
           hjust = 0) +
  scale_x_continuous(breaks = seq(0, 100, 25)) +
  scale_y_continuous(breaks = seq(0, 100, 25)) +
  theme(text = element_text(size = 20))

# ggsave(str_c("../../figures/plots/exp3_", model, "_scatter.pdf"),
#        width = 5,
#        height = 5)

4.5.2.4 Ternary plots for individual regression results

df.plot = df.exp3.causal.individual_fit %>% 
  select(participant, fit_whs) %>% 
  mutate(estimates = map(fit_whs, tidy)) %>% 
  select(participant, estimates) %>% 
  unnest(estimates) %>% 
  filter(str_detect(term, "b_"),
         !str_detect(term, "Intercept")) %>% 
  mutate(term = str_remove(term, "b_")) %>% 
  select(participant, term, estimate) %>% 
  spread(term, estimate) %>% 
  mutate_at(vars(how:whether), .funs = list(norm = ~ . / (how + whether + sufficient))) %>% 
  mutate(color = 0,
         color = ifelse(how_norm == max(how_norm), 1, color),
         color = ifelse(whether_norm == max(whether_norm), 2, color),
         color = ifelse(sufficient_norm == max(sufficient_norm), 3, color),
         color = factor(color))

ggplot(data = df.plot,
       mapping = aes(x = how,
                     y = sufficient,
                     z = whether,
                     color = color)) + 
  geom_point(alpha = 0.7,
             size = 2,
             show.legend = F) + 
  scale_color_manual(values = c("black", "red", "green", "blue")) + 
  coord_tern() +
  theme_bw() + 
  theme_showarrows() + 
  theme(text = element_text(size = 20),
        tern.axis.title.T = element_blank(),
        tern.axis.title.L = element_blank(),
        tern.axis.title.R = element_blank())

# ggsave(str_c("../../figures/plots/exp3_individual_regression_ternary_plot_scaled.pdf"),
#        width = 5,
#        height = 5)

4.5.2.5 Individual participant variance for a selection of clips (clustered)

# \u2713 = check mark 
# \u2717 = cross mark 

check = "\u2713"
cross = "\u2717"

x_labels = c(str_c("\nW ", cross, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross))

df.plot = df.exp3.causal.long %>%
  filter(clip %in% c(7, 23, 3, 15, 16)) %>%
  mutate(clip = factor(clip, levels = c(7, 23, 3, 15, 16),
    labels = c("causal chain", "double prevention", "joint causation",
      "overdetermination", "preemption")))

df.cluster = df.exp3.causal.individual_fit %>%
  select(participant, fit_whs) %>%
  mutate(estimates = map(fit_whs, tidy)) %>%
  select(participant, estimates) %>%
  unnest(estimates) %>%
  filter(str_detect(term, "b_"),
         !str_detect(term, "Intercept")) %>%
  mutate(term = str_remove(term, "b_")) %>%
  select(participant, term, estimate) %>%
  spread(term, estimate) %>%
  mutate(group = ifelse(how > whether, "how", "whether"))

df.plot = df.plot %>% 
  left_join(df.cluster %>% 
              select(participant, group),
            by = "participant")

df.labels = df.plot %>% 
  distinct(clip, ball) %>% 
  arrange(clip, ball) %>% 
  mutate(labels = x_labels,
         rating = -10,
         group = NA,
         participant = NA)

func_load_image = function(clip){
  readPNG(str_c("../../figures/diagrams/exp3_clip", clip, ".png"))
}

df.clips = df.plot %>% 
  distinct(clip) %>% 
  arrange(clip) %>% 
  mutate(number = c(7, 23, 3, 15, 16),
         grob = map(.x = number, .f = ~ func_load_image(clip = .x)),
         group = NA,
         participant = NA,
         ball = NA,
         label = str_c("clip ", number))

ball_a = readPNG("../../figures/diagrams/ballA.png")
ball_a = rasterGrob(ball_a, interpolate = TRUE)

ball_b = readPNG("../../figures/diagrams/ballB.png")
ball_b = rasterGrob(ball_b, interpolate = TRUE)

# ggplot(df.plot,aes(x=ball,y=rating,group=id,color = best))+
p = ggplot(data = df.plot, 
       mapping = aes(x = ball,
                     y = rating,
                     group = participant,
                     shape = group)) +
  geom_line(mapping = aes(linetype = "individual",
                          color = group),
            size = 1,
            alpha = 0.3) +
  geom_point(mapping = aes(color = group),
             size = 1,
             alpha = 0.3) +
  stat_summary(mapping = aes(group = group,
                             color = group,
                             linetype = "mean"),
               fun.y = "mean",
               geom = "line",
               size = 1.5) +
  stat_summary(data = df.plot %>% filter(group == "how"),
               mapping = aes(group = group,
                             color = group),
               fun.data = "mean_cl_boot",
               geom = "pointrange",
               size = 1,
               shape = 19) +
  stat_summary(data = df.plot %>% filter(group == "whether"),
               mapping = aes(group = group,
                             color = group),
               fun.data = "mean_cl_boot",
               geom = "pointrange",
               size = 1,
               shape = 19) +
  annotation_custom(grob = ball_a, xmin = 0.5, xmax = 1.5, ymin = -30, ymax = -10) + 
  annotation_custom(grob = ball_b, xmin = 1.5, xmax = 2.5, ymin = -30, ymax = -10) +
  geom_text(data = df.labels,
            mapping = aes(label = labels,
                          y = -Inf),
            vjust = 1.2,
            family = "Arial Unicode MS",
            size = 8) + 
  geom_custom(data = df.clips,
              mapping = aes(data = grob, x = 1.5, y = Inf),
              grob_fun = function(x) rasterGrob(x,
                                                interpolate = T,
                                                vjust = -0.25)) +
  geom_text(data = df.clips,
            mapping = aes(label = label,
                          y = Inf,
                          x = -Inf),
            size = 9,
            hjust = -0.2,
            vjust = -4) + 
  facet_wrap(~clip,
             ncol = 8) +
  coord_cartesian(xlim = c(0.9, 2.1),
                  ylim = c(-5, 105),
                  expand = F,
                  clip = "off") +
  scale_y_continuous(breaks = seq(0, 100, 25),
                     labels = seq(0, 100, 25)) +
  scale_color_brewer(palette = "Set1",
                     guide = "none") +
  labs(y = "causal responsibility rating",
       linetype = "legend",
       shape = "") +
  theme_bw() +
  guides(linetype = guide_legend(override.aes = list(alpha = c(0.3, 1)),
                                 keywidth = unit(1.2, "cm")),
         shape = guide_legend(override.aes = list(color = c("#e41a1c", "#377eb8"),
                                                  shape = c(19, 19),
                                                  alpha = c(1, 1),
                                                  size = c(5, 5)))) +
  theme(panel.grid = element_blank(),
        text = element_text(size = 20),
        legend.position = c(0.505, 0.23),
        axis.title.x = element_blank(),
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 24),
        legend.background = element_rect(fill = "transparent"),
        strip.text = element_text(size = 30),
        axis.text.x = element_blank(),
        legend.key = element_rect(fill = "transparent"),
        legend.box = "horizontal",
        legend.spacing.x = unit(0.1, "cm"),
        panel.spacing.x = unit(1, "cm"),
        plot.margin = margin(b = 5, l = 0.2, r = 0.2, t = 7.5, unit = "cm"))

# ggsave(str_c("../../figures/plots/exp3_individual_variance_selection_lines_clustered.pdf"),
#        plot = p,
#        width = 20,
#        height = 8.5,
#        device = cairo_pdf)

4.5.3 Tables

4.5.3.1 Relationship between predictors

df.exp3.model %>% 
  filter(noise == unique(df.exp3.model.counterfactual$noise)) %>%
  gather("index", "value", A_difference:B_robust) %>% 
  separate(index, into = c("ball", "aspect")) %>% 
  spread(aspect, value) %>% 
  select(difference, 
         whether, 
         how,
         sufficient,
         robust) %>% 
  correlate() %>% 
  shave() %>% 
  fashion() %>% 
  xtable() %>% 
  print(include.rownames = F,
        booktabs = T)

Correlation method: 'pearson'
Missing treated using: 'pairwise.complete.obs'

% latex table generated in R 3.6.1 by xtable 1.8-4 package
% Tue Nov 12 15:06:10 2019
\begin{table}[ht]
\centering
\begin{tabular}{llllll}
  \toprule
rowname & difference & whether & how & sufficient & robust \\ 
  \midrule
difference &  &  &  &  &  \\ 
  whether &  .48 &  &  &  &  \\ 
  how &  .76 &  .24 &  &  &  \\ 
  sufficient &  .18 &  .05 &  .34 &  &  \\ 
  robust &  .42 &  .94 &  .21 &  .13 &  \\ 
   \bottomrule
\end{tabular}
\end{table}

4.5.3.2 Population level predictors in the mixed effects models

fit.brm_exp3_causal_w %>% 
  tidy() %>% 
  filter(str_detect(term, "b_")) %>% 
  mutate(model = "CSM_w") %>% 
  bind_rows(fit.brm_exp3_causal_wh %>% 
              tidy() %>% 
              filter(str_detect(term, "b_")) %>% 
              mutate(model = "CSM_wh")) %>% 
  bind_rows(fit.brm_exp3_causal_whs %>% 
              tidy() %>% 
              filter(str_detect(term, "b_")) %>% 
              mutate(model = "CSM_whs")) %>% 
  mutate(term = tolower(term),
         term = str_remove_all(term, "b_")) %>% 
  rename(`lower 95% HDI` = lower,
         `upper 95% HDI` = upper) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  select(model, everything()) %>% 
  xtable() %>% 
  print(include.rownames = F,
        booktabs = T)

% latex table generated in R 3.6.1 by xtable 1.8-4 package
% Tue Nov 12 15:06:11 2019
\begin{table}[ht]
\centering
\begin{tabular}{llrrrr}
  \toprule
model & term & estimate & std.error & lower 95\% HDI & upper 95\% HDI \\ 
  \midrule
CSM\_w & intercept & 31.68 & 1.61 & 29.02 & 34.37 \\ 
  CSM\_w & whether & 39.30 & 2.16 & 35.76 & 42.82 \\ 
  CSM\_wh & intercept & 10.29 & 1.56 & 7.80 & 12.83 \\ 
  CSM\_wh & whether & 28.40 & 2.68 & 23.96 & 32.70 \\ 
  CSM\_wh & how & 36.07 & 2.60 & 31.80 & 40.27 \\ 
  CSM\_whs & intercept & 10.43 & 1.59 & 7.82 & 13.01 \\ 
  CSM\_whs & whether & 27.93 & 2.60 & 23.67 & 32.18 \\ 
  CSM\_whs & how & 26.47 & 2.91 & 21.76 & 31.35 \\ 
  CSM\_whs & sufficient & 30.90 & 3.06 & 25.77 & 35.96 \\ 
   \bottomrule
\end{tabular}
\end{table}

4.5.3.3 Individual participants regression results

df.exp3.causal.individual_fit %>% 
  select_if(~ (is.numeric(.) | is.factor(.))) %>% 
  mutate_at(.vars = vars(contains("_w")), .funs = ~ round(., 2)) %>% 
  select(participant, everything(), best_model) %>% 
  xtable() %>% 
  print(include.rownames = F,
        booktabs = T)

% latex table generated in R 3.6.1 by xtable 1.8-4 package
% Tue Nov 12 15:06:12 2019
\begin{table}[ht]
\centering
\begin{tabular}{rrrrrrrl}
  \toprule
participant & r\_w & r\_wh & r\_whs & rmse\_w & rmse\_wh & rmse\_whs & best\_model \\ 
  \midrule
  1 & 0.29 & 0.38 & 0.48 & 30.00 & 29.03 & 27.74 & whs \\ 
    2 & 0.20 & 0.77 & 0.77 & 39.32 & 27.94 & 27.55 & whs \\ 
    3 & 0.37 & 0.78 & 0.79 & 38.90 & 28.00 & 26.95 & whs \\ 
    4 & 0.41 & 0.55 & 0.63 & 43.62 & 40.60 & 38.31 & whs \\ 
    5 & 0.45 & 0.53 & 0.54 & 39.23 & 37.32 & 36.81 & whs \\ 
    6 & 0.18 & 0.54 & 0.54 & 40.84 & 36.32 & 36.07 & whs \\ 
    7 & 0.49 & 0.61 & 0.64 & 32.96 & 29.93 & 29.15 & whs \\ 
    8 & 0.39 & 0.63 & 0.68 & 38.37 & 33.26 & 31.64 & whs \\ 
    9 & 0.40 & 0.72 & 0.76 & 35.86 & 28.34 & 26.39 & whs \\ 
   10 & 0.28 & 0.52 & 0.53 & 29.13 & 26.28 & 25.85 & whs \\ 
   11 & 0.51 & 0.56 & 0.60 & 33.66 & 32.19 & 31.17 & whs \\ 
   12 & 0.56 & 0.57 & 0.70 & 32.78 & 32.08 & 28.84 & whs \\ 
   13 & 0.56 & 0.71 & 0.74 & 29.83 & 25.40 & 24.21 & whs \\ 
   14 & 0.43 & 0.47 & 0.50 & 33.72 & 32.82 & 32.23 & whs \\ 
   15 & 0.19 & 0.35 & 0.37 & 39.77 & 38.29 & 37.94 & whs \\ 
   16 & 0.62 & 0.67 & 0.76 & 40.49 & 37.98 & 34.39 & whs \\ 
   17 & 0.42 & 0.66 & 0.71 & 28.57 & 23.82 & 22.42 & whs \\ 
   18 & 0.31 & 0.35 & 0.39 & 37.78 & 37.09 & 36.57 & whs \\ 
   19 & 0.37 & 0.60 & 0.62 & 33.41 & 29.23 & 28.62 & whs \\ 
   20 & 0.57 & 0.61 & 0.68 & 36.47 & 35.00 & 32.71 & whs \\ 
   21 & 0.43 & 0.50 & 0.57 & 30.95 & 29.61 & 28.32 & whs \\ 
   22 & 0.51 & 0.67 & 0.71 & 33.72 & 29.36 & 28.00 & whs \\ 
   23 & 0.33 & 0.45 & 0.52 & 38.76 & 37.00 & 35.57 & whs \\ 
   24 & 0.46 & 0.63 & 0.69 & 34.01 & 29.96 & 28.22 & whs \\ 
   25 & 0.61 & 0.66 & 0.70 & 30.60 & 28.91 & 27.55 & whs \\ 
   26 & 0.57 & 0.71 & 0.76 & 40.43 & 35.03 & 32.80 & whs \\ 
   27 & 0.46 & 0.52 & 0.66 & 43.12 & 41.49 & 38.12 & whs \\ 
   28 & 0.30 & 0.45 & 0.46 & 29.33 & 27.62 & 27.38 & whs \\ 
   29 & 0.48 & 0.60 & 0.65 & 38.93 & 36.03 & 34.48 & whs \\ 
   30 & 0.18 & 0.70 & 0.70 & 33.06 & 24.95 & 24.71 & whs \\ 
   31 & 0.48 & 0.59 & 0.69 & 37.54 & 34.72 & 31.76 & whs \\ 
   32 & 0.37 & 0.61 & 0.70 & 40.62 & 35.94 & 32.92 & whs \\ 
   33 & 0.50 & 0.56 & 0.64 & 38.64 & 36.79 & 34.85 & whs \\ 
   34 & 0.24 & 0.51 & 0.52 & 44.50 & 40.86 & 40.19 & whs \\ 
   35 & 0.15 & 0.32 & 0.30 & 34.56 & 33.36 & 33.37 & wh \\ 
   36 & 0.19 & 0.84 & 0.82 & 45.38 & 29.58 & 29.70 & wh \\ 
   37 & 0.55 & 0.63 & 0.67 & 36.04 & 33.36 & 32.11 & whs \\ 
   38 & 0.52 & 0.61 & 0.67 & 36.03 & 33.39 & 31.53 & whs \\ 
   39 & 0.47 & 0.77 & 0.77 & 29.61 & 21.75 & 21.58 & whs \\ 
   40 & 0.46 & 0.78 & 0.79 & 32.41 & 23.60 & 22.69 & whs \\ 
   41 & 0.32 & 0.50 & 0.58 & 32.90 & 30.38 & 28.69 & whs \\ 
   \bottomrule
\end{tabular}
\end{table}

4.5.3.4 CSMwhs predictions for selection of cases

df.table = df.exp3.model.causal %>%
  left_join(df.exp3.causal.regression) %>% 
  filter(clip %in% c(7, 23, 16, 3, 15)) %>%
  mutate_at(.vars = vars(difference, robust, sufficient, whether, heuristic), .funs = ~ round(., 2)) %>%
  select(clip, ball, difference, whether, how, sufficient, robust) %>%
  mutate(clip = factor(clip, levels = c(7, 23, 16, 3, 15))) %>%
  mutate_at(.vars = vars(difference, whether, how, sufficient, robust), .funs = ~ ifelse(. < 0.5, paste0("xmark (", ., ")"), paste0("cmark (", ., ")"))) %>%
  arrange(clip)

Joining, by = c("clip", "ball")

xtable(df.table) %>%
  print(include.rownames = FALSE)

% latex table generated in R 3.6.1 by xtable 1.8-4 package
% Tue Nov 12 15:06:12 2019
\begin{table}[ht]
\centering
\begin{tabular}{lllllll}
  \hline
clip & ball & difference & whether & how & sufficient & robust \\ 
  \hline
7 & A & cmark (1) & xmark (0.2) & cmark (1) & xmark (0) & xmark (0.18) \\ 
  7 & B & cmark (1) & cmark (1) & cmark (1) & cmark (0.81) & cmark (0.6) \\ 
  23 & A & xmark (0.05) & xmark (0) & xmark (0) & xmark (0) & xmark (0) \\ 
  23 & B & cmark (0.96) & cmark (0.89) & xmark (0) & xmark (0) & cmark (0.8) \\ 
  16 & A & cmark (1) & xmark (0.19) & cmark (1) & cmark (1) & xmark (0.29) \\ 
  16 & B & xmark (0) & xmark (0) & xmark (0) & xmark (0) & xmark (0) \\ 
  3 & A & cmark (1) & cmark (0.88) & cmark (1) & xmark (0.13) & cmark (0.74) \\ 
  3 & B & cmark (1) & cmark (0.87) & cmark (1) & xmark (0.12) & cmark (0.74) \\ 
  15 & A & cmark (1) & xmark (0) & cmark (1) & cmark (1) & xmark (0.08) \\ 
  15 & B & cmark (1) & xmark (0) & cmark (1) & cmark (1) & xmark (0.06) \\ 
   \hline
\end{tabular}
\end{table}

4.5.3.5 CSMwhs predictions for all cases

df.table = df.exp3.causal.regression %>% 
  left_join(df.exp3.model.causal,
            by = c("clip", "ball")) %>% 
  left_join(df.exp3.clipinfo %>% 
              select(-clipindex),
            by = c("clip")) %>% 
  mutate_at(.vars = vars(difference, whether, how, sufficient, robust),
            .funs = ~ . * 100) %>% 
  mutate_at(.vars = vars(-ball), .funs = ~ round(., 0)) %>%
  select(clip, ball, contains("outcome"), difference, whether, how, sufficient, robust, w, wh, whs, heuristic, rating = rating_mean)

xtable(df.table, digits = 0) %>%
  print(include.rownames = FALSE)

% latex table generated in R 3.6.1 by xtable 1.8-4 package
% Tue Nov 12 15:06:12 2019
\begin{table}[ht]
\centering
\begin{tabular}{rlrrrrrrrrrrrrrr}
  \hline
clip & ball & outcome\_both & outcome\_a & outcome\_b & outcome\_none & difference & whether & how & sufficient & robust & w & wh & whs & heuristic & rating \\ 
  \hline
1 & A & 0 & 0 & 0 & 0 & 100 & 41 & 100 & 20 & 40 & 48 & 58 & 54 & 75 & 42 \\ 
  1 & B & 0 & 0 & 0 & 0 & 100 & 19 & 100 & 16 & 7 & 39 & 52 & 47 & 72 & 37 \\ 
  2 & A & 0 & 0 & 0 & 0 & 58 & 11 & 0 & 0 & 12 & 36 & 13 & 14 & 64 & 21 \\ 
  2 & B & 0 & 0 & 0 & 0 & 14 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 64 & 19 \\ 
  3 & A & 1 & 0 & 0 & 0 & 100 & 88 & 100 & 13 & 74 & 66 & 71 & 66 & 92 & 76 \\ 
  3 & B & 1 & 0 & 0 & 0 & 100 & 87 & 100 & 12 & 74 & 66 & 71 & 65 & 92 & 75 \\ 
  4 & A & 1 & 0 & 0 & 0 & 100 & 76 & 100 & 0 & 70 & 62 & 68 & 58 & 76 & 63 \\ 
  4 & B & 1 & 0 & 0 & 0 & 100 & 99 & 100 & 16 & 60 & 71 & 74 & 69 & 72 & 78 \\ 
  5 & A & 0 & 0 & 1 & 0 & 100 & 97 & 100 & 0 & 49 & 70 & 74 & 64 & 68 & 71 \\ 
  5 & B & 0 & 0 & 1 & 0 & 100 & 0 & 100 & 0 & 0 & 32 & 46 & 37 & 88 & 22 \\ 
  6 & A & 0 & 0 & 1 & 0 & 100 & 67 & 100 & 12 & 38 & 58 & 66 & 59 & 71 & 73 \\ 
  6 & B & 0 & 0 & 1 & 0 & 100 & 12 & 100 & 4 & 12 & 36 & 50 & 42 & 72 & 22 \\ 
  7 & A & 1 & 0 & 1 & 0 & 100 & 20 & 100 & 0 & 18 & 40 & 52 & 43 & 86 & 59 \\ 
  7 & B & 1 & 0 & 1 & 0 & 100 & 100 & 100 & 81 & 60 & 71 & 75 & 90 & 85 & 79 \\ 
  8 & A & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 65 & 7 \\ 
  8 & B & 1 & 0 & 1 & 0 & 100 & 100 & 100 & 100 & 100 & 71 & 75 & 96 & 98 & 92 \\ 
  9 & A & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 53 & 8 \\ 
  9 & B & 0 & 1 & 0 & 0 & 100 & 100 & 100 & 0 & 100 & 71 & 75 & 65 & 90 & 90 \\ 
  10 & A & 0 & 1 & 0 & 0 & 82 & 13 & 0 & 0 & 17 & 37 & 14 & 14 & 55 & 23 \\ 
  10 & B & 0 & 1 & 0 & 0 & 100 & 88 & 0 & 0 & 68 & 66 & 35 & 35 & 59 & 55 \\ 
  11 & A & 1 & 1 & 0 & 0 & 100 & 74 & 100 & 76 & 70 & 61 & 67 & 81 & 82 & 93 \\ 
  11 & B & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 54 & 4 \\ 
  12 & A & 1 & 1 & 0 & 0 & 100 & 84 & 100 & 75 & 81 & 65 & 70 & 84 & 75 & 77 \\ 
  12 & B & 1 & 1 & 0 & 0 & 100 & 0 & 100 & 9 & 25 & 32 & 46 & 40 & 72 & 37 \\ 
  13 & A & 0 & 1 & 1 & 0 & 71 & 44 & 0 & 0 & 40 & 49 & 23 & 23 & 54 & 8 \\ 
  13 & B & 0 & 1 & 1 & 0 & 81 & 51 & 0 & 0 & 47 & 52 & 25 & 25 & 59 & 64 \\ 
  14 & A & 0 & 1 & 1 & 0 & 100 & 97 & 0 & 0 & 60 & 70 & 38 & 38 & 59 & 22 \\ 
  14 & B & 0 & 1 & 1 & 0 & 96 & 86 & 0 & 0 & 58 & 66 & 35 & 34 & 58 & 18 \\ 
  15 & A & 1 & 1 & 1 & 0 & 100 & 0 & 100 & 100 & 8 & 32 & 46 & 68 & 91 & 76 \\ 
  15 & B & 1 & 1 & 1 & 0 & 100 & 0 & 100 & 100 & 6 & 32 & 46 & 68 & 91 & 76 \\ 
  16 & A & 1 & 1 & 1 & 0 & 100 & 19 & 100 & 100 & 29 & 39 & 52 & 73 & 94 & 92 \\ 
  16 & B & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 53 & 4 \\ 
  17 & A & 0 & 0 & 0 & 1 & 100 & 11 & 100 & 53 & 14 & 36 & 49 & 56 & 86 & 69 \\ 
  17 & B & 0 & 0 & 0 & 1 & 100 & 0 & 100 & 45 & 18 & 32 & 46 & 51 & 84 & 46 \\ 
  18 & A & 0 & 0 & 0 & 1 & 100 & 6 & 100 & 52 & 12 & 34 & 48 & 54 & 74 & 63 \\ 
  18 & B & 0 & 0 & 0 & 1 & 100 & 4 & 100 & 50 & 6 & 33 & 48 & 54 & 74 & 66 \\ 
  19 & A & 1 & 0 & 0 & 1 & 100 & 78 & 100 & 3 & 71 & 62 & 69 & 60 & 73 & 53 \\ 
  19 & B & 1 & 0 & 0 & 1 & 100 & 80 & 100 & 4 & 71 & 63 & 69 & 60 & 73 & 49 \\ 
  20 & A & 1 & 0 & 0 & 1 & 100 & 92 & 100 & 3 & 72 & 68 & 73 & 64 & 75 & 41 \\ 
  20 & B & 1 & 0 & 0 & 1 & 100 & 94 & 100 & 3 & 59 & 69 & 73 & 64 & 74 & 71 \\ 
  21 & A & 0 & 0 & 1 & 1 & 100 & 56 & 100 & 53 & 53 & 54 & 62 & 69 & 76 & 80 \\ 
  21 & B & 0 & 0 & 1 & 1 & 100 & 4 & 100 & 24 & 4 & 33 & 47 & 45 & 78 & 18 \\ 
  22 & A & 0 & 0 & 1 & 1 & 100 & 100 & 100 & 96 & 84 & 71 & 75 & 95 & 69 & 60 \\ 
  22 & B & 0 & 0 & 1 & 1 & 100 & 4 & 100 & 0 & 14 & 33 & 48 & 38 & 71 & 42 \\ 
  23 & A & 1 & 0 & 1 & 1 & 5 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 54 & 3 \\ 
  23 & B & 1 & 0 & 1 & 1 & 96 & 89 & 0 & 0 & 80 & 67 & 35 & 35 & 60 & 39 \\ 
  24 & A & 1 & 0 & 1 & 1 & 100 & 61 & 100 & 1 & 54 & 56 & 64 & 54 & 75 & 44 \\ 
  24 & B & 1 & 0 & 1 & 1 & 100 & 98 & 100 & 21 & 84 & 70 & 74 & 71 & 72 & 73 \\ 
  25 & A & 0 & 1 & 0 & 1 & 100 & 18 & 100 & 17 & 19 & 39 & 52 & 47 & 80 & 43 \\ 
  25 & B & 0 & 1 & 0 & 1 & 100 & 78 & 100 & 64 & 69 & 62 & 68 & 78 & 75 & 73 \\ 
  26 & A & 0 & 1 & 0 & 1 & 100 & 10 & 100 & 3 & 12 & 36 & 49 & 41 & 78 & 39 \\ 
  26 & B & 0 & 1 & 0 & 1 & 100 & 96 & 100 & 42 & 57 & 69 & 74 & 77 & 68 & 69 \\ 
  27 & A & 1 & 1 & 0 & 1 & 100 & 100 & 100 & 47 & 98 & 71 & 75 & 79 & 78 & 80 \\ 
  27 & B & 1 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 55 & 6 \\ 
  28 & A & 1 & 1 & 0 & 1 & 100 & 94 & 100 & 14 & 84 & 69 & 73 & 68 & 92 & 89 \\ 
  28 & B & 1 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 53 & 5 \\ 
  29 & A & 0 & 1 & 1 & 1 & 100 & 64 & 100 & 21 & 49 & 57 & 65 & 61 & 78 & 47 \\ 
  29 & B & 0 & 1 & 1 & 1 & 100 & 70 & 100 & 23 & 43 & 59 & 66 & 64 & 73 & 67 \\ 
  30 & A & 0 & 1 & 1 & 1 & 100 & 65 & 100 & 28 & 55 & 57 & 65 & 64 & 83 & 58 \\ 
  30 & B & 0 & 1 & 1 & 1 & 100 & 55 & 100 & 24 & 45 & 53 & 62 & 60 & 78 & 56 \\ 
  31 & A & 1 & 1 & 1 & 1 & 100 & 0 & 100 & 1 & 0 & 32 & 46 & 37 & 75 & 44 \\ 
  31 & B & 1 & 1 & 1 & 1 & 100 & 0 & 100 & 1 & 0 & 32 & 47 & 37 & 71 & 46 \\ 
  32 & A & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 54 & 5 \\ 
  32 & B & 1 & 1 & 1 & 1 & 100 & 66 & 100 & 52 & 65 & 58 & 65 & 71 & 77 & 71 \\ 
   \hline
\end{tabular}
\end{table}

4.5.3.6 Heuristic model

fit.brm_exp3_causal_heuristic %>% 
  tidy() %>% 
  filter(str_detect(term, "b_")) %>% 
  mutate(term = str_remove(term, "b_"),
         term = tolower(term)) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  rename(`lower 95% HDI` = lower,
         `upper 95% HDI` = upper) %>% 
  xtable() %>%
  print(include.rownames = F,
        booktabs = T)

% latex table generated in R 3.6.1 by xtable 1.8-4 package
% Tue Nov 12 15:06:12 2019
\begin{table}[ht]
\centering
\begin{tabular}{lrrrr}
  \toprule
term & estimate & std.error & lower 95\% HDI & upper 95\% HDI \\ 
  \midrule
intercept & 49.79 & 1.46 & 47.39 & 52.23 \\ 
  moving & 0.22 & 0.22 & 0.01 & 0.66 \\ 
  speed & 2.06 & 0.85 & 0.64 & 3.43 \\ 
  contact\_e & 0.38 & 0.35 & 0.02 & 1.07 \\ 
  e\_speed\_diff & 0.12 & 0.12 & 0.01 & 0.37 \\ 
  e\_direction\_diff & 1.04 & 0.74 & 0.08 & 2.43 \\ 
  total\_speed\_diff & 2.21 & 0.95 & 0.67 & 3.80 \\ 
  total\_direction\_diff & 3.99 & 0.91 & 2.45 & 5.45 \\ 
  transfer & 15.59 & 0.79 & 14.30 & 16.90 \\ 
  e\_moving & 0.18 & 0.17 & 0.01 & 0.51 \\ 
  exclusive & 4.39 & 0.73 & 3.16 & 5.56 \\ 
   \bottomrule
\end{tabular}
\end{table}

5 Session info

sessionInfo()

R version 3.6.1 (2019-07-05)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.6

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] forcats_0.4.0     stringr_1.4.0     dplyr_0.8.3      
 [4] purrr_0.3.2       readr_1.3.1       tidyr_1.0.0      
 [7] tibble_2.1.3      tidyverse_1.2.1   egg_0.4.5        
[10] gridExtra_2.3     png_0.1-7         tidybayes_1.1.0  
[13] Hmisc_4.2-0       Formula_1.2-3     survival_2.44-1.1
[16] lattice_0.20-38   ggtern_3.1.0      ggplot2_3.2.1    
[19] janitor_1.2.0     jsonlite_1.6      xtable_1.8-4     
[22] corrr_0.4.0       brms_2.10.0       Rcpp_1.0.2       
[25] broom_0.5.2       lme4_1.1-21       Matrix_1.2-17    
[28] knitr_1.25       

loaded via a namespace (and not attached):
  [1] readxl_1.3.1              backports_1.1.4          
  [3] plyr_1.8.4                igraph_1.2.4.1           
  [5] lazyeval_0.2.2            splines_3.6.1            
  [7] svUnit_0.7-12             crosstalk_1.0.0          
  [9] rstantools_2.0.0          inline_0.3.15            
 [11] digest_0.6.21             htmltools_0.3.6          
 [13] rsconnect_0.8.15          fansi_0.4.0              
 [15] magrittr_1.5              checkmate_1.9.4          
 [17] cluster_2.1.0             modelr_0.1.5             
 [19] bayesm_3.1-3              matrixStats_0.55.0       
 [21] xts_0.11-2                prettyunits_1.0.2        
 [23] colorspace_1.4-1          rvest_0.3.4              
 [25] haven_2.1.1               xfun_0.9                 
 [27] callr_3.3.2               crayon_1.3.4             
 [29] zeallot_0.1.0             zoo_1.8-6                
 [31] glue_1.3.1                gtable_0.3.0             
 [33] compositions_1.40-2       pkgbuild_1.0.5           
 [35] rstan_2.19.2              DEoptimR_1.0-8           
 [37] abind_1.4-5               scales_1.0.0             
 [39] miniUI_0.1.1.1            htmlTable_1.13.2         
 [41] HDInterval_0.2.0          ggstance_0.3.3           
 [43] foreign_0.8-72            latex2exp_0.4.0          
 [45] stats4_3.6.1              StanHeaders_2.19.0       
 [47] DT_0.9                    htmlwidgets_1.3          
 [49] httr_1.4.1                threejs_0.3.1            
 [51] arrayhelpers_1.0-20160527 RColorBrewer_1.1-2       
 [53] ellipsis_0.3.0            acepack_1.4.1            
 [55] pkgconfig_2.0.3           loo_2.1.0                
 [57] nnet_7.3-12               utf8_1.1.4               
 [59] labeling_0.3              tidyselect_0.2.5         
 [61] rlang_0.4.0               reshape2_1.4.3           
 [63] later_0.8.0               cellranger_1.1.0         
 [65] munsell_0.5.0             tools_3.6.1              
 [67] cli_1.1.0                 generics_0.0.2           
 [69] ggridges_0.5.1            evaluate_0.14            
 [71] yaml_2.2.0                processx_3.4.1           
 [73] robustbase_0.93-5         nlme_3.1-141             
 [75] mime_0.7                  xml2_1.2.2               
 [77] compiler_3.6.1            bayesplot_1.7.0          
 [79] shinythemes_1.1.2         rstudioapi_0.10          
 [81] stringi_1.4.3             ps_1.3.0                 
 [83] Brobdingnag_1.2-6         nloptr_1.2.1             
 [85] markdown_1.1              shinyjs_1.0              
 [87] tensorA_0.36.1            vctrs_0.2.0              
 [89] pillar_1.4.2              lifecycle_0.1.0          
 [91] bridgesampling_0.7-2      data.table_1.12.2        
 [93] httpuv_1.5.2              R6_2.4.0                 
 [95] latticeExtra_0.6-28       bookdown_0.13            
 [97] promises_1.0.1            boot_1.3-23              
 [99] energy_1.7-6              colourpicker_1.0         
[101] MASS_7.3-51.4             gtools_3.8.1             
[103] assertthat_0.2.1          proto_1.0.0              
[105] withr_2.1.2               shinystan_2.5.0          
[107] parallel_3.6.1            hms_0.5.1                
[109] rpart_4.1-15              coda_0.19-3              
[111] minqa_1.2.4               snakecase_0.11.0         
[113] rmarkdown_1.15            shiny_1.3.2              
[115] lubridate_1.7.4           base64enc_0.1-3          
[117] dygraphs_1.1.1.6